From 927f0c20673081aa600776231f7c3dafff4b5eb1 Mon Sep 17 00:00:00 2001
From: Johanna <johanna@zijderveld.de>
Date: Mon, 18 Dec 2023 15:17:15 +0100
Subject: [PATCH] hermitian cost function graphene example (untested and
 non-generic dimensional)

---
 examples/graphene_extended_hubbard.ipynb | 262 +++++++++++++++++++----
 1 file changed, 224 insertions(+), 38 deletions(-)

diff --git a/examples/graphene_extended_hubbard.ipynb b/examples/graphene_extended_hubbard.ipynb
index e5aa8f6..c645c71 100644
--- a/examples/graphene_extended_hubbard.ipynb
+++ b/examples/graphene_extended_hubbard.ipynb
@@ -7,7 +7,16 @@
    "metadata": {
     "tags": []
    },
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n"
+     ]
+    }
+   ],
    "source": [
     "import kwant\n",
     "import numpy as np\n",
@@ -30,9 +39,75 @@
   {
    "cell_type": "code",
    "execution_count": 2,
-   "id": "9ebc3cce-0338-4616-8021-9fecee76f178",
    "metadata": {},
    "outputs": [],
+   "source": [
+    "def flat_to_matrix(flat, shape):\n",
+    "    matrix = np.zeros(shape, dtype=complex)\n",
+    "    matrix[:, :, *np.triu_indices(shape[-1])] = flat.reshape(shape[0], shape[1], -1)\n",
+    "    lower_triangle_wo_diag = matrix.transpose(0, 1, 3, 2)[\n",
+    "        :, :, *np.tril_indices(shape[-1], k=-1)\n",
+    "    ]\n",
+    "    matrix[:, :, *np.tril_indices(shape[-1], k=-1)] = lower_triangle_wo_diag.reshape(\n",
+    "        shape[0], shape[1], -1\n",
+    "    ).conj()\n",
+    "\n",
+    "    return matrix\n",
+    "\n",
+    "def matrix_to_flat(matrix):\n",
+    "    matrix[\n",
+    "        :,\n",
+    "        :,\n",
+    "        *np.triu_indices(matrix.shape[-1]),\n",
+    "    ].flatten()\n",
+    "\n",
+    "def upper_triangle_cost(delta_mf_flatten, model):\n",
+    "    mfk_shape = model.mf_k.shape\n",
+    "\n",
+    "    # From flat to non-flat:\n",
+    "    mf_delta_shaped = flat_to_matrix(delta_mf_flatten, mfk_shape)\n",
+    "\n",
+    "    # Doing the update\n",
+    "    mf_shaped = model.mf_k + mf_delta_shaped\n",
+    "    _, model.mf_k = hf.updated_matrices(mf_k=mf_shaped, model=model)\n",
+    "    delta_mf = mf_shaped - model.mf_k\n",
+    "\n",
+    "    # Now flattening again\n",
+    "    return matrix_to_flat(delta_mf)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "9ebc3cce-0338-4616-8021-9fecee76f178",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
+      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n"
+     ]
+    }
+   ],
    "source": [
     "# Create translationally-invariant `kwant.Builder`\n",
     "graphene_builder, int_builder = kwant_examples.graphene_extended_hubbard()"
@@ -48,7 +123,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "id": "c2dd6c3c-d9bb-4833-b1a0-d28c5d5a935c",
    "metadata": {},
    "outputs": [],
@@ -66,46 +141,59 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "params = dict(U=1, V=1)\n",
+    "int_model = utils.builder2tb_model(int_builder, params)\n",
+    "model = hf.Model(tb_model, int_model=int_model)\n",
+    "vectors = utils.generate_vectors(10, model.dim)\n",
+    "model.vectors = [*vectors, *model.tb_model.keys()]\n",
+    "model.random_guess(model.vectors)\n",
+    "model.kgrid_evaluation(nk=15)\n",
+    "matrix_to_flat(model.mf_k)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
    "id": "41bd9f60-8f29-4e7c-a0c4-a0bbf66445b2",
    "metadata": {},
    "outputs": [],
    "source": [
     "def compute_gap(\n",
-    "    tb_model,\n",
-    "    int_model,\n",
+    "    model,\n",
     "    nk,\n",
     "    nk_dense,\n",
-    "    U,\n",
     "    filling=2,\n",
     "    guess=None,\n",
     "):\n",
-    "    scale = np.max(np.array([*tb_model.values()]))\n",
-    "    # Generate guess on the same grid\n",
-    "    if guess is None:\n",
-    "        guess = utils.generate_guess(nk, tb_model, int_model, scale=scale)\n",
-    "    else:\n",
-    "        guess += utils.generate_guess(\n",
-    "            nk, tb_model, int_model, scale=scale\n",
-    "        )\n",
+    "\n",
+    "    model.guess = guess\n",
     "\n",
     "    # Find groundstate Hamiltonian on the same grid\n",
-    "    mf_model, mf = hf.find_groundstate_ham(\n",
-    "        tb_model=tb_model,\n",
-    "        int_model=int_model,\n",
+    "    mf_model = hf.find_groundstate_ham(\n",
+    "        model,\n",
     "        filling=filling,\n",
     "        nk=nk,\n",
-    "        guess=guess,\n",
-    "        return_mf=True,\n",
+    "        cutoff_Vk=10,  # no intuition about cutoff_Vk\n",
+    "        cost_function=upper_triangle_cost,\n",
+    "    )\n",
+    "\n",
+    "    mf_k = utils.kgrid_hamiltonian(  # only used for diagonalization to get gap\n",
+    "        nk=nk_dense,\n",
+    "        hk=utils.model2hk(tb_model=mf_model),\n",
+    "        dim=2\n",
     "    )\n",
-    "    # Compute groundstate Hamiltonian on a dense grid\n",
-    "    scf_ham = utils.kgrid_hamiltonian(nk_dense, mf_model)\n",
     "    # Diagonalize groundstate Hamiltonian\n",
-    "    vals, vecs = np.linalg.eigh(scf_ham)\n",
+    "    vals, _ = np.linalg.eigh(mf_k)\n",
     "    # Extract dense-grid Fermi energy\n",
     "    E_F = utils.get_fermi_energy(vals, filling)\n",
     "    gap = utils.calc_gap(vals, E_F)\n",
-    "    return gap, mf"
+    "\n",
+    "    # the guess was kind of unclear\n",
+    "    return gap, mf_model"
    ]
   },
   {
@@ -118,7 +206,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 6,
    "id": "32b9e7c5-db12-44f9-930c-21e5494404b8",
    "metadata": {
     "tags": []
@@ -132,12 +220,13 @@
     "        gap_U = []\n",
     "        for V in Vs:\n",
     "            params = dict(U=U, V=V)\n",
+    "            int_model = utils.builder2tb_model(int_builder, params)\n",
+    "            model = hf.Model(tb_model, int_model=int_model)\n",
     "            _gap, guess = compute_gap(\n",
-    "                tb_model=tb_model,\n",
-    "                int_model=utils.builder2tb_model(int_builder, params),\n",
+    "                model=model,\n",
     "                nk=nk,\n",
     "                nk_dense=nk_dense,\n",
-    "                guess=guess, U=U\n",
+    "                guess=guess,\n",
     "            )\n",
     "            gap_U.append(_gap)\n",
     "        gap.append(gap_U)\n",
@@ -154,7 +243,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 14,
    "id": "6a8c08a9-7e31-420b-b6b4-709abfb26793",
    "metadata": {
     "tags": []
@@ -164,28 +253,95 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 20/20 [53:16<00:00, 159.83s/it]\n"
+      "  0%|          | 0/10 [00:00<?, ?it/s]\n"
+     ]
+    },
+    {
+     "ename": "ValueError",
+     "evalue": "cannot reshape array of size 1 into shape (15,15,newaxis)",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
+      "\u001b[1;32m/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb Cell 13\u001b[0m line \u001b[0;36m7\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=4'>5</a>\u001b[0m Us \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mlinspace(\u001b[39m0\u001b[39m, \u001b[39m3\u001b[39m, \u001b[39m10\u001b[39m, endpoint\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m)\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=5'>6</a>\u001b[0m Vs \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mlinspace(\u001b[39m0\u001b[39m, \u001b[39m1.5\u001b[39m, \u001b[39m10\u001b[39m, endpoint\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m)\n\u001b[0;32m----> <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=6'>7</a>\u001b[0m gap \u001b[39m=\u001b[39m compute_phase_diagram(Us, Vs, nk\u001b[39m=\u001b[39;49mnk, nk_dense\u001b[39m=\u001b[39;49mnk_dense)\n",
+      "\u001b[1;32m/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb Cell 13\u001b[0m line \u001b[0;36m1\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=7'>8</a>\u001b[0m     int_model \u001b[39m=\u001b[39m utils\u001b[39m.\u001b[39mbuilder2tb_model(int_builder, params)\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=8'>9</a>\u001b[0m     model \u001b[39m=\u001b[39m hf\u001b[39m.\u001b[39mModel(tb_model, int_model\u001b[39m=\u001b[39mint_model)\n\u001b[0;32m---> <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=9'>10</a>\u001b[0m     _gap, guess \u001b[39m=\u001b[39m compute_gap(\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=10'>11</a>\u001b[0m         model\u001b[39m=\u001b[39;49mmodel,\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=11'>12</a>\u001b[0m         nk\u001b[39m=\u001b[39;49mnk,\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=12'>13</a>\u001b[0m         nk_dense\u001b[39m=\u001b[39;49mnk_dense,\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=13'>14</a>\u001b[0m         guess\u001b[39m=\u001b[39;49mguess,\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=14'>15</a>\u001b[0m     )\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=15'>16</a>\u001b[0m     gap_U\u001b[39m.\u001b[39mappend(_gap)\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=16'>17</a>\u001b[0m gap\u001b[39m.\u001b[39mappend(gap_U)\n",
+      "\u001b[1;32m/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb Cell 13\u001b[0m line \u001b[0;36m1\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=8'>9</a>\u001b[0m model\u001b[39m.\u001b[39mguess \u001b[39m=\u001b[39m guess\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=10'>11</a>\u001b[0m \u001b[39m# Find groundstate Hamiltonian on the same grid\u001b[39;00m\n\u001b[0;32m---> <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=11'>12</a>\u001b[0m mf_model \u001b[39m=\u001b[39m hf\u001b[39m.\u001b[39;49mfind_groundstate_ham(\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=12'>13</a>\u001b[0m     model,\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=13'>14</a>\u001b[0m     filling\u001b[39m=\u001b[39;49mfilling,\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=14'>15</a>\u001b[0m     nk\u001b[39m=\u001b[39;49mnk,\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=15'>16</a>\u001b[0m     cutoff_Vk\u001b[39m=\u001b[39;49m\u001b[39m10\u001b[39;49m,  \u001b[39m# no intuition about cutoff_Vk\u001b[39;49;00m\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=16'>17</a>\u001b[0m     cost_function\u001b[39m=\u001b[39;49mupper_triangle_cost,\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=17'>18</a>\u001b[0m )\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=19'>20</a>\u001b[0m mf_k \u001b[39m=\u001b[39m utils\u001b[39m.\u001b[39mkgrid_hamiltonian(  \u001b[39m# only used for diagonalization to get gap\u001b[39;00m\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=20'>21</a>\u001b[0m     nk\u001b[39m=\u001b[39mnk_dense,\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=21'>22</a>\u001b[0m     hk\u001b[39m=\u001b[39mutils\u001b[39m.\u001b[39mmodel2hk(tb_model\u001b[39m=\u001b[39mmf_model),\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=22'>23</a>\u001b[0m     dim\u001b[39m=\u001b[39m\u001b[39m2\u001b[39m\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=23'>24</a>\u001b[0m )\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=24'>25</a>\u001b[0m \u001b[39m# Diagonalize groundstate Hamiltonian\u001b[39;00m\n",
+      "File \u001b[0;32m~/kwant-scf/examples/codes/hf.py:250\u001b[0m, in \u001b[0;36mfind_groundstate_ham\u001b[0;34m(model, cutoff_Vk, filling, nk, cost_function, optimizer, optimizer_kwargs)\u001b[0m\n\u001b[1;32m    248\u001b[0m model\u001b[39m.\u001b[39mkgrid_evaluation(nk\u001b[39m=\u001b[39mnk)\n\u001b[1;32m    249\u001b[0m fun \u001b[39m=\u001b[39m partial(cost_function, model\u001b[39m=\u001b[39mmodel)\n\u001b[0;32m--> 250\u001b[0m mf_k \u001b[39m=\u001b[39m optimizer(fun, matrix_to_flat(model\u001b[39m.\u001b[39;49mmf_k), \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49moptimizer_kwargs)\n\u001b[1;32m    251\u001b[0m \u001b[39massert\u001b[39;00m np\u001b[39m.\u001b[39mallclose((mf_k \u001b[39m-\u001b[39m np\u001b[39m.\u001b[39mmoveaxis(mf_k, \u001b[39m-\u001b[39m\u001b[39m1\u001b[39m, \u001b[39m-\u001b[39m\u001b[39m2\u001b[39m)\u001b[39m.\u001b[39mconj()) \u001b[39m/\u001b[39m \u001b[39m2\u001b[39m, \u001b[39m0\u001b[39m, atol\u001b[39m=\u001b[39m\u001b[39m1e-5\u001b[39m)\n\u001b[1;32m    252\u001b[0m mf_k \u001b[39m=\u001b[39m (mf_k \u001b[39m+\u001b[39m np\u001b[39m.\u001b[39mmoveaxis(mf_k, \u001b[39m-\u001b[39m\u001b[39m1\u001b[39m, \u001b[39m-\u001b[39m\u001b[39m2\u001b[39m)\u001b[39m.\u001b[39mconj()) \u001b[39m/\u001b[39m \u001b[39m2\u001b[39m\n",
+      "File \u001b[0;32m<string>:6\u001b[0m, in \u001b[0;36manderson\u001b[0;34m(F, xin, iter, alpha, w0, M, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, **kw)\u001b[0m\n",
+      "File \u001b[0;32m~/opt/anaconda3/envs/python3/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:172\u001b[0m, in \u001b[0;36mnonlin_solve\u001b[0;34m(F, x0, jacobian, iter, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, full_output, raise_exception)\u001b[0m\n\u001b[1;32m    169\u001b[0m x \u001b[39m=\u001b[39m x0\u001b[39m.\u001b[39mflatten()\n\u001b[1;32m    171\u001b[0m dx \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mfull_like(x, np\u001b[39m.\u001b[39minf)\n\u001b[0;32m--> 172\u001b[0m Fx \u001b[39m=\u001b[39m func(x)\n\u001b[1;32m    173\u001b[0m Fx_norm \u001b[39m=\u001b[39m norm(Fx)\n\u001b[1;32m    175\u001b[0m jacobian \u001b[39m=\u001b[39m asjacobian(jacobian)\n",
+      "File \u001b[0;32m~/opt/anaconda3/envs/python3/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:168\u001b[0m, in \u001b[0;36mnonlin_solve.<locals>.func\u001b[0;34m(z)\u001b[0m\n\u001b[1;32m    167\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mfunc\u001b[39m(z):\n\u001b[0;32m--> 168\u001b[0m     \u001b[39mreturn\u001b[39;00m _as_inexact(F(_array_like(z, x0)))\u001b[39m.\u001b[39mflatten()\n",
+      "\u001b[1;32m/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb Cell 13\u001b[0m line \u001b[0;36m2\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=20'>21</a>\u001b[0m mfk_shape \u001b[39m=\u001b[39m model\u001b[39m.\u001b[39mmf_k\u001b[39m.\u001b[39mshape\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=22'>23</a>\u001b[0m \u001b[39m# From flat to non-flat:\u001b[39;00m\n\u001b[0;32m---> <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=23'>24</a>\u001b[0m mf_delta_shaped \u001b[39m=\u001b[39m flat_to_matrix(delta_mf_flatten, mfk_shape)\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=25'>26</a>\u001b[0m \u001b[39m# Doing the update\u001b[39;00m\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=26'>27</a>\u001b[0m mf_shaped \u001b[39m=\u001b[39m model\u001b[39m.\u001b[39mmf_k \u001b[39m+\u001b[39m mf_delta_shaped\n",
+      "\u001b[1;32m/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb Cell 13\u001b[0m line \u001b[0;36m3\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mflat_to_matrix\u001b[39m(flat, shape):\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=1'>2</a>\u001b[0m     matrix \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mzeros(shape, dtype\u001b[39m=\u001b[39m\u001b[39mcomplex\u001b[39m)\n\u001b[0;32m----> <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m     matrix[:, :, \u001b[39m*\u001b[39mnp\u001b[39m.\u001b[39mtriu_indices(shape[\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m])] \u001b[39m=\u001b[39m flat\u001b[39m.\u001b[39;49mreshape(shape[\u001b[39m0\u001b[39;49m], shape[\u001b[39m1\u001b[39;49m], \u001b[39m-\u001b[39;49m\u001b[39m1\u001b[39;49m)\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m     lower_triangle_wo_diag \u001b[39m=\u001b[39m matrix\u001b[39m.\u001b[39mtranspose(\u001b[39m0\u001b[39m, \u001b[39m1\u001b[39m, \u001b[39m3\u001b[39m, \u001b[39m2\u001b[39m)[\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=4'>5</a>\u001b[0m         :, :, \u001b[39m*\u001b[39mnp\u001b[39m.\u001b[39mtril_indices(shape[\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m], k\u001b[39m=\u001b[39m\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m)\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=5'>6</a>\u001b[0m     ]\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=6'>7</a>\u001b[0m     matrix[:, :, \u001b[39m*\u001b[39mnp\u001b[39m.\u001b[39mtril_indices(shape[\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m], k\u001b[39m=\u001b[39m\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m)] \u001b[39m=\u001b[39m lower_triangle_wo_diag\u001b[39m.\u001b[39mreshape(\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=7'>8</a>\u001b[0m         shape[\u001b[39m0\u001b[39m], shape[\u001b[39m1\u001b[39m], \u001b[39m-\u001b[39m\u001b[39m1\u001b[39m\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/rzijderveld/kwant-scf/examples/graphene_extended_hubbard.ipynb#X13sZmlsZQ%3D%3D?line=8'>9</a>\u001b[0m     )\u001b[39m.\u001b[39mconj()\n",
+      "\u001b[0;31mValueError\u001b[0m: cannot reshape array of size 1 into shape (15,15,newaxis)"
      ]
     }
    ],
    "source": [
     "# Generate dense-grid k-points\n",
     "# Interaction strengths\n",
-    "nk=15\n",
-    "nk_dense=30\n",
-    "Us = np.linspace(0, 3, 20, endpoint=True)\n",
-    "Vs = np.linspace(0, 1.5, 20, endpoint=True)\n",
+    "nk = 15\n",
+    "nk_dense = 30\n",
+    "Us = np.linspace(0, 3, 10, endpoint=True)\n",
+    "Vs = np.linspace(0, 1.5, 10, endpoint=True)\n",
     "gap = compute_phase_diagram(Us, Vs, nk=nk, nk_dense=nk_dense)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "> \u001b[0;32m/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_80069/2781680018.py\u001b[0m(3)\u001b[0;36mflat_to_matrix\u001b[0;34m()\u001b[0m\n",
+      "\u001b[0;32m      1 \u001b[0;31m\u001b[0;32mdef\u001b[0m \u001b[0mflat_to_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mflat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m      2 \u001b[0;31m    \u001b[0mmatrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcomplex\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m----> 3 \u001b[0;31m    \u001b[0mmatrix\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtriu_indices\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mflat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m      4 \u001b[0;31m    lower_triangle_wo_diag = matrix.transpose(0, 1, 3, 2)[\n",
+      "\u001b[0m\u001b[0;32m      5 \u001b[0;31m        \u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtril_indices\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\n",
+      "> \u001b[0;32m/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_80069/2781680018.py\u001b[0m(24)\u001b[0;36mupper_triangle_cost\u001b[0;34m()\u001b[0m\n",
+      "\u001b[0;32m     22 \u001b[0;31m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m     23 \u001b[0;31m    \u001b[0;31m# From flat to non-flat:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m---> 24 \u001b[0;31m    \u001b[0mmf_delta_shaped\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mflat_to_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdelta_mf_flatten\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmfk_shape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m     25 \u001b[0;31m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m     26 \u001b[0;31m    \u001b[0;31m# Doing the update\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\n",
+      "array(nan)\n",
+      "> \u001b[0;32m/Users/rzijderveld/opt/anaconda3/envs/python3/lib/python3.11/site-packages/scipy/optimize/_nonlin.py\u001b[0m(168)\u001b[0;36mfunc\u001b[0;34m()\u001b[0m\n",
+      "\u001b[0;32m    166 \u001b[0;31m    \u001b[0mx0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_as_inexact\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m    167 \u001b[0;31m    \u001b[0;32mdef\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m--> 168 \u001b[0;31m        \u001b[0;32mreturn\u001b[0m \u001b[0m_as_inexact\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mF\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_array_like\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m    169 \u001b[0;31m    \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m    170 \u001b[0;31m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\n",
+      "array([nan])\n",
+      "array(nan)\n",
+      "> \u001b[0;32m/Users/rzijderveld/opt/anaconda3/envs/python3/lib/python3.11/site-packages/scipy/optimize/_nonlin.py\u001b[0m(172)\u001b[0;36mnonlin_solve\u001b[0;34m()\u001b[0m\n",
+      "\u001b[0;32m    170 \u001b[0;31m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m    171 \u001b[0;31m    \u001b[0mdx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfull_like\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m--> 172 \u001b[0;31m    \u001b[0mFx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m    173 \u001b[0;31m    \u001b[0mFx_norm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mFx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\u001b[0;32m    174 \u001b[0;31m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0m\n"
+     ]
+    }
+   ],
+   "source": [
+    "%debug"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 45,
    "id": "e17fc96c-c463-4e1f-8250-c254d761b92a",
    "metadata": {},
    "outputs": [],
    "source": [
     "import xarray as xr\n",
+    "\n",
     "gap_da = xr.DataArray(data=gap, coords=dict(Us=Us, Vs=Vs))"
    ]
   },
@@ -197,6 +353,36 @@
     "We note that the gap openings coincide with the phase transitions from gapless to charge density wave or antiferromagnetic groundstates as predicted in [arXiv:1204.4531](https://arxiv.org/abs/1204.4531). "
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 46,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.collections.QuadMesh at 0x14e254450>"
+      ]
+     },
+     "execution_count": 46,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGwCAYAAAAJ/wd3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5q0lEQVR4nO3df1Rd1Z3//9cFAqRRsAmGJOaGUDUJyphJL+0IkfojCS50suyvCWNmQlSYhiExQ6j6CWY+NeZjpVZLsTVgMiYyWRNTliZqXaWau76aHwb9LqGkdYytVlPBeAkDbSGJBuTe8/kjcj/ecEEu98c5l/t8zNprhs3e7Pfuacd3997nbJthGIYAAAAiJM7sAAAAQGwh+QAAABFF8gEAACKK5AMAAEQUyQcAAIgokg8AABBRJB8AACCiEswOINI8Ho8++ugjXXjhhbLZbGaHAwCwMMMwdOrUKc2aNUtxceH77+tnz57VwMBA0H8nMTFRycnJIYgovGIu+fjoo49kt9vNDgMAEEU6Ojo0e/bssPzts2fPKjPjAnV2uYP+WzNmzNDx48ctn4DEXPJx4YUXSpLy/u4eJSQkmRwNAMDKBgf71fz//9j7z45wGBgYUGeXWx+0zlXKheNfXek75VGG408aGBgg+bCaoa2WhIQkJSRY++EAAKwhEtv0F1xo0wUXjn8cj6LnKAEHTgEAsAC34Qm6jEddXZ0yMzOVnJwsh8Ohw4cPj9p+9+7dWrhwob70pS9p5syZuv3229XT0xPQmCQfAABYgEdG0CVQjY2Nqqio0KZNm9TW1qb8/HwVFhaqvb3db/tXX31VxcXFKikp0VtvvaWnn35ab7zxhkpLSwMal+QDAIAYVVNTo5KSEpWWliorK0u1tbWy2+2qr6/32/7111/X3LlztX79emVmZuqaa67RmjVr1NLSEtC4JB8AAFiAJwT/I0l9fX0+pb+/3+94AwMDam1tVUFBgU99QUGBmpub/fbJy8vThx9+qKamJhmGoZMnT+qZZ57RzTffHNBcST4AALAAt2EEXSTJbrcrNTXVW6qrq/2O193dLbfbrfT0dJ/69PR0dXZ2+u2Tl5en3bt3q6ioSImJiZoxY4Yuuugi/fznPw9oriQfAABMIB0dHert7fWWqqqqUduf/yaPYRgjvt1z7NgxrV+/Xj/4wQ/U2tqqF198UcePH1dZWVlAMcbcq7YAAFjReA+Nfr6/JKWkpCglJeUL26elpSk+Pn7YKkdXV9ew1ZAh1dXVWrx4se6++25J0lVXXaUpU6YoPz9fDzzwgGbOnDmmWFn5AADAAjwy5A6iBJq4JCYmyuFwyOl0+tQ7nU7l5eX57fPxxx8P+8x8fHy8pHMrJmNF8gEAQIyqrKzUE088oZ07d+rtt9/Whg0b1N7e7t1GqaqqUnFxsbf98uXLtW/fPtXX1+v999/XkSNHtH79en3961/XrFmzxjwu2y4AAFhAqLZdAlFUVKSenh5t2bJFLpdL2dnZampqUkZGhiTJ5XL5fPPjtttu06lTp/TYY4/p+9//vi666CLdcMMNeuihhwIa12YEsk4yAfT19Sk1NVXfWPy/+bw6AGBUg4NndejI/1Fvb++YzlGMx9A/l955O10XBnG3y6lTHs3LOhnWWEOFbRcAABBRbLsAAGABns9KMP2jBckHAAAWMPTWSjD9owXJBwAAFuA2zpVg+kcLznwAAICIYuUDAAAL4MwHAACIKI9scsv/nSpj7R8t2HYBAAARxcoHAAAW4DHOlWD6RwuSDwAALMAd5LZLMH0jjW0XAAAQUax8AABgAbG08kHyAQCABXgMmzxGEG+7BNE30th2AQAAEcXKBwAAFsC2CwAAiCi34uQOYkPCHcJYwo3kAwAACzCCPPNhcOYDAADAP1Y+AACwAM58AACAiHIbcXIbQZz5iKLPq7PtAgAAIoqVDwAALMAjmzxBrAl4FD1LHyQfAABYAGc+IuTQoUN6+OGH1draKpfLpWeffVbf/OY3x9T3yJEjuvbaa5Wdna2jR48GPPakP3+shPhoeisaABBpNne/2SFMSKae+Thz5owWLlyoxx57LKB+vb29Ki4u1pIlS8IUGQAAkTV04DSYEi1MXfkoLCxUYWFhwP3WrFmjlStXKj4+Xs8991zoAwMAIMLOnfkI4mK5KNp2iZ406TNPPvmk3nvvPd13331jat/f36++vj6fAgAAzBNVB07fffddbdy4UYcPH1ZCwthCr66u1v333x/myAAACI4nyLtdoultl6hZ+XC73Vq5cqXuv/9+zZs3b8z9qqqq1Nvb6y0dHR1hjBIAgPHhzIcFnTp1Si0tLWpra9O6deskSR6PR4ZhKCEhQfv379cNN9wwrF9SUpKSkpIiHS4AAAHxKI7vfFhNSkqK3nzzTZ+6uro6vfzyy3rmmWeUmZlpUmQAACAQpiYfp0+f1h//+Efvz8ePH9fRo0c1depUzZkzR1VVVTpx4oR27dqluLg4ZWdn+/SfPn26kpOTh9UDABBt3IZNbiOIj4wF0TfSTE0+WlpadP3113t/rqyslCStXr1aDQ0Ncrlcam9vNys8AAAixh3kgVM32y5jc91118kwRv4Xq6GhYdT+mzdv1ubNm0MbFAAACKvoORoLAMAE5jHigi7jUVdXp8zMTCUnJ8vhcOjw4cMjtr3ttttks9mGlSuvvDKgMUk+AACwgKFtl2BKoBobG1VRUaFNmzapra1N+fn5KiwsHPHIw6OPPiqXy+UtHR0dmjp1qv7hH/4hoHFJPgAAiFE1NTUqKSlRaWmpsrKyVFtbK7vdrvr6er/tU1NTNWPGDG9paWnRX/7yF91+++0BjRs1r9oCADCReRTcGyuez/73+deIjPS9q4GBAbW2tmrjxo0+9QUFBWpubh7TmDt27NDSpUuVkZERUKysfAAAYAFDHxkLpkiS3W5Xamqqt1RXV/sdr7u7W263W+np6T716enp6uzs/MJ4XS6Xfv3rX6u0tDTgubLyAQDABNLR0aGUlBTvz1/0lW+bzXe1xTCMYXX+NDQ06KKLLtI3v/nNgGMk+QAAwAKCvZ9lqG9KSopP8jGStLQ0xcfHD1vl6OrqGrYacj7DMLRz506tWrVKiYmJAcfKtgsAABbgkS3oEojExEQ5HA45nU6feqfTqby8vFH7Hjx4UH/84x9VUlIS8DwlVj4AALCEUK18BKKyslKrVq1STk6OcnNztX37drW3t6usrEySfK45+bwdO3bo7/7u78Z9vQnJBwAAMaqoqEg9PT3asmWLXC6XsrOz1dTU5H17xd81J729vdq7d68effTRcY9L8gEAgAUEf7fL+PqWl5ervLzc7+/8XXOSmpqqjz/+eFxjDSH5AADAAjyGTZ5gvvMRRbfacuAUAABEFCsfAABYgCfIbRdPFK0nkHwAAGABwdxMO9Q/WkRPpAAAYEJg5QMAAAtwyyZ3gB8KO79/tCD5AADAAth2AQAACBNWPgAAsAC3gts6cYculLAj+QAAwAJiaduF5AMAAAsw42I5s0RPpAAAYEJg5QMAAAswZJMniDMfBq/aAgCAQLDtAgAAECYxu/Lxl4VTFZ+YbHYYAAALcw+cld6OzFgewyaPMf6tk2D6RlrMJh8AAFiJO8hbbYPpG2nREykAAJgQWPkAAMAC2HYBAAAR5VGcPEFsSATTN9KiJ1IAADAhsPIBAIAFuA2b3EFsnQTTN9JIPgAAsIBYOvNh6rbLoUOHtHz5cs2aNUs2m03PPffcqO337dunZcuW6eKLL1ZKSopyc3P10ksvRSZYAADCyPjsVtvxFoMvnI7NmTNntHDhQj322GNjan/o0CEtW7ZMTU1Nam1t1fXXX6/ly5erra0tzJECAIBQMXXbpbCwUIWFhWNuX1tb6/Pzgw8+qOeff14vvPCCFi1a5LdPf3+/+vv7vT/39fWNK1YAAMLJLZvcQVwOF0zfSIueNRo/PB6PTp06palTp47Yprq6Wqmpqd5it9sjGCEAAGPjMf7fuY/xFbNnMHZRnXz85Cc/0ZkzZ7RixYoR21RVVam3t9dbOjo6IhghAAA4X9S+7bJnzx5t3rxZzz//vKZPnz5iu6SkJCUlJUUwMgAAAjd0cDSY/tEiKpOPxsZGlZSU6Omnn9bSpUvNDgcAgKB5ZJMniHMbwfSNtOhJkz6zZ88e3XbbbXrqqad08803mx0OAAAIkKkrH6dPn9Yf//hH78/Hjx/X0aNHNXXqVM2ZM0dVVVU6ceKEdu3aJelc4lFcXKxHH31UV199tTo7OyVJkydPVmpqqilzAAAgFGLpC6emrny0tLRo0aJF3tdkKysrtWjRIv3gBz+QJLlcLrW3t3vbb9u2TYODg1q7dq1mzpzpLf/2b/9mSvwAAIRKMB8YC/a8SKSZuvJx3XXXyTBGfjeooaHB5+cDBw6ENyAAABB2UXngFACAicajIO924cApAAAIhPHZ2y7jLcY4k4+6ujplZmYqOTlZDodDhw8fHrV9f3+/Nm3apIyMDCUlJenSSy/Vzp07AxqTlQ8AACzAjFttGxsbVVFRobq6Oi1evFjbtm1TYWGhjh07pjlz5vjts2LFCp08eVI7duzQZZddpq6uLg0ODgY0LskHAAATyPl3mI32sc2amhqVlJSotLRU0rk71F566SXV19erurp6WPsXX3xRBw8e1Pvvv++92mTu3LkBx8i2CwAAFhCqt13sdrvPnWb+kghJGhgYUGtrqwoKCnzqCwoK1Nzc7LfPL3/5S+Xk5OjHP/6xLrnkEs2bN0933XWXPvnkk4DmysoHAAAWEKptl46ODqWkpHjrR1r16O7ultvtVnp6uk99enq69zta53v//ff16quvKjk5Wc8++6y6u7tVXl6uP//5zwGd+yD5AABgAklJSfFJPr6Izeab8BiGMaxuiMfjkc1m0+7du70f96ypqdF3v/tdbd26VZMnTx7TmGy7AABgAcG86TKee2HS0tIUHx8/bJWjq6tr2GrIkJkzZ+qSSy7x+ap4VlaWDMPQhx9+OOaxST4AALCAoW2XYEogEhMT5XA45HQ6feqdTqfy8vL89lm8eLE++ugjnT592lv3zjvvKC4uTrNnzx7z2CQfAADEqMrKSj3xxBPauXOn3n77bW3YsEHt7e0qKyuTJFVVVam4uNjbfuXKlZo2bZpuv/12HTt2TIcOHdLdd9+tO+64Y8xbLhJnPgAAsAQzvvNRVFSknp4ebdmyRS6XS9nZ2WpqalJGRoak4XesXXDBBXI6nbrzzjuVk5OjadOmacWKFXrggQcCGpfkAwAACzAj+ZCk8vJylZeX+/3d+XesSdKCBQuGbdUEim0XAAAQUax8AABgAWatfJiB5AMAAAswFNzNtEboQgk7kg8AACwgllY+OPMBAAAiipUPAAAsIJZWPmI2+YjvN5TgiaYdMgAYv+6/sfZCd9qbHrND8O/TyP1zIpaSD2v/uxEAAEw4MbvyAQCAlcTSygfJBwAAFmAYNhlBJBDB9I00tl0AAEBEsfIBAIAFeGQL6iNjwfSNNJIPAAAsIJbOfLDtAgAAIoqVDwAALCCWDpySfAAAYAGxtO1C8gEAgAXE0soHZz4AAEBEsfIBAIAFGEFuu0TTygfJBwAAFmBIMoK4xy6arkpl2wUAAESUqcnHoUOHtHz5cs2aNUs2m03PPffcF/Y5ePCgHA6HkpOT9ZWvfEWPP/54+AMFACDMhr5wGkyJFqYmH2fOnNHChQv12GOPjan98ePHddNNNyk/P19tbW269957tX79eu3duzfMkQIAEF5Db7sEU6KFqWc+CgsLVVhYOOb2jz/+uObMmaPa2lpJUlZWllpaWvTII4/oO9/5jt8+/f396u/v9/7c19cXVMwAACA4UXXm47XXXlNBQYFP3Y033qiWlhZ9+umnfvtUV1crNTXVW+x2eyRCBQAgIEMfGQumRIuoSj46OzuVnp7uU5eenq7BwUF1d3f77VNVVaXe3l5v6ejoiESoAAAExDCCL9Ei6l61tdl8Mzvjs3+1z68fkpSUpKSkpLDHBQAAxiaqko8ZM2aos7PTp66rq0sJCQmaNm2aSVEBABC8WPq8elQlH7m5uXrhhRd86vbv36+cnBxNmjTJpKgAAAheLCUfpp75OH36tI4ePaqjR49KOvcq7dGjR9Xe3i7p3HmN4uJib/uysjJ98MEHqqys1Ntvv62dO3dqx44duuuuu8wIHwCAkImlA6emrny0tLTo+uuv9/5cWVkpSVq9erUaGhrkcrm8iYgkZWZmqqmpSRs2bNDWrVs1a9Ys/exnPxvxNVsAAGA9piYf1113nffAqD8NDQ3D6q699lr95je/CWNUAABEXrBvrPC2CwAACMi55COYMx8hDCbMouo7HwAAIPqRfAAAYAFm3e1SV1enzMxMJScny+Fw6PDhwyO2PXDggGw227Dy+9//PqAx2XYBAMACjM9KMP0D1djYqIqKCtXV1Wnx4sXatm2bCgsLdezYMc2ZM2fEfn/4wx+UkpLi/fniiy8OaFxWPgAAiFE1NTUqKSlRaWmpsrKyVFtbK7vdrvr6+lH7TZ8+XTNmzPCW+Pj4gMYl+QAAwAJCte3S19fnUz5/s/vnDQwMqLW1ddiFrQUFBWpubh411kWLFmnmzJlasmSJXnnllYDnSvIBAIAVGCEokux2u89t7tXV1X6H6+7ultvt9nth6/lXmQyZOXOmtm/frr1792rfvn2aP3++lixZokOHDgU0Vc58AABgBUF+Xl2f9e3o6PA5j/FFl6v6u7B1pMta58+fr/nz53t/zs3NVUdHhx555BF94xvfGHOorHwAADCBpKSk+JSRko+0tDTFx8f7vbD1/NWQ0Vx99dV69913A4qR5AMAAAsY+sJpMCUQiYmJcjgccjqdPvVOp1N5eXlj/jttbW2aOXNmQGOz7QIAgAWYcattZWWlVq1apZycHOXm5mr79u1qb29XWVmZpHMXvJ44cUK7du2SJNXW1mru3Lm68sorNTAwoP/6r//S3r17tXfv3oDGjdnkI/X1DiXEJZodBgBExEfXjfzNBitIPdL+xY1MMOgZMDuEsCoqKlJPT4+2bNkil8ul7OxsNTU1KSMjQ5KGXfA6MDCgu+66SydOnNDkyZN15ZVX6le/+pVuuummgMaN2eQDAABLMWzeQ6Pj7j8O5eXlKi8v9/u78y94veeee3TPPfeMa5zPI/kAAMACYulWWw6cAgCAiGLlAwAAKzDjcheTkHwAAGABZrztYha2XQAAQESx8gEAgFVE0dZJMEg+AACwgFjadiH5AADACmLowClnPgAAQESx8gEAgCXYPivB9I8OJB8AAFgB2y4AAADhwcoHAABWEEMrHyQfAABYgUm32pqBbRcAABBRrHwAAGABhnGuBNM/WpB8AABgBTF05oNtFwAAEFGsfAAAYAUxdOCU5AMAAAuwGedKMP2jhenbLnV1dcrMzFRycrIcDocOHz48avvdu3dr4cKF+tKXvqSZM2fq9ttvV09PT4SiBQAgTIwQlChhavLR2NioiooKbdq0SW1tbcrPz1dhYaHa29v9tn/11VdVXFyskpISvfXWW3r66af1xhtvqLS0NMKRAwCA8TI1+aipqVFJSYlKS0uVlZWl2tpa2e121dfX+23/+uuva+7cuVq/fr0yMzN1zTXXaM2aNWppaYlw5AAAhNjQmY9gSpQwLfkYGBhQa2urCgoKfOoLCgrU3Nzst09eXp4+/PBDNTU1yTAMnTx5Us8884xuvvnmEcfp7+9XX1+fTwEAwHLYdgm/7u5uud1upaen+9Snp6ers7PTb5+8vDzt3r1bRUVFSkxM1IwZM3TRRRfp5z//+YjjVFdXKzU11VvsdntI5wEAAAJj+oFTm813mcgwjGF1Q44dO6b169frBz/4gVpbW/Xiiy/q+PHjKisrG/HvV1VVqbe311s6OjpCGj8AACERQysfpr1qm5aWpvj4+GGrHF1dXcNWQ4ZUV1dr8eLFuvvuuyVJV111laZMmaL8/Hw98MADmjlz5rA+SUlJSkpKCv0EAAAIJb5wGn6JiYlyOBxyOp0+9U6nU3l5eX77fPzxx4qL8w05Pj5e0rkVEwAAYH2mfmSssrJSq1atUk5OjnJzc7V9+3a1t7d7t1Gqqqp04sQJ7dq1S5K0fPly/cu//Ivq6+t14403yuVyqaKiQl//+tc1a9YsM6cCAEBw+MJpZBQVFamnp0dbtmyRy+VSdna2mpqalJGRIUlyuVw+3/y47bbbdOrUKT322GP6/ve/r4suukg33HCDHnroIbOmAABASMTSF04DTj7+8z//U2lpad7XW++55x5t375dV1xxhfbs2eNNHMaqvLxc5eXlfn/X0NAwrO7OO+/UnXfeGWjYAADAIgI+8/Hggw9q8uTJkqTXXntNjz32mH784x8rLS1NGzZsCHmAAADEBN52GVlHR4cuu+wySdJzzz2n7373u/re976nxYsX67rrrgt1fAAAYIIJeOXjggsu8F7ktn//fi1dulSSlJycrE8++SS00QEAECNs+n/nPsZVzJ5AAAJe+Vi2bJlKS0u1aNEivfPOO96zH2+99Zbmzp0b6vgAAMAEM+bk4+jRo/rbv/1bbd26Vf/+7/+ujo4O7d27V9OmTZMktba26tZbbw1boKE28JXp8iQkmx0GAETEtLl/MTuEUQ1cNsPsEPwaHDwruSI0mEmv2tbV1enhhx+Wy+XSlVdeqdraWuXn539hvyNHjujaa69Vdna2jh49GtCYY04+vvrVr2rRokUqLS3VD3/4Q6Wmpvr8/v777w9oYAAA8DkmfOG0sbFRFRUVqqur0+LFi7Vt2zYVFhbq2LFjmjNnzoj9ent7VVxcrCVLlujkyZMBjzvmMx9HjhzRV7/6VW3cuFEzZ87UqlWr9MorrwQ8IAAACJ/zb3Lv7+8fsW1NTY1KSkpUWlqqrKws1dbWym63q76+ftQx1qxZo5UrVyo3N3dcMY45+cjNzdV//Md/qLOzU/X19ero6NDSpUt16aWX6oc//KE+/PDDcQUAAAAUsldt7Xa7z23u1dXVfocbGBhQa2urCgoKfOoLCgrU3Nw8YphPPvmk3nvvPd13333jnmrAb7tMnjxZq1ev1oEDB/TOO+/o1ltv1bZt25SZmambbrpp3IEAABDLgnrT5XNfR+3o6PC5zb2qqsrveN3d3XK73cMuc01PTx926euQd999Vxs3btTu3buVkDD+j6QH9Xn1Sy+9VBs3bpTdbte9996rl156KZg/BwAAgpSSkqKUlJQxt7fZfA+qGoYxrE6S3G63Vq5cqfvvv1/z5s0LKsZxJx8HDx7Uzp07tXfvXsXHx2vFihUqKSkJKhgAAGJWhA+cpqWlKT4+ftgqR1dX17DVEEk6deqUWlpa1NbWpnXr1kmSPB6PDMNQQkKC9u/frxtuuGFMYweUfHR0dKihoUENDQ06fvy48vLy9POf/1wrVqzQlClTAvlTAADg8yKcfCQmJsrhcMjpdOpb3/qWt97pdOqWW24Z1j4lJUVvvvmmT11dXZ1efvllPfPMM8rMzBzz2GNOPpYtW6ZXXnlFF198sYqLi3XHHXdo/vz5Yx4IAABYS2VlpVatWqWcnBzl5uZq+/btam9vV1lZmSSpqqpKJ06c0K5duxQXF6fs7Gyf/tOnT1dycvKw+i8y5uRj8uTJ2rt3r/7+7/9e8fHxAQ0CAABG9/lDo+PtH6iioiL19PRoy5Ytcrlcys7OVlNTk/eGepfLpfb29vEHNQKbYRhRdA9e8Pr6+pSamqpvLP7fSuALpwBixF//1xmzQxjVRQ9Zc+t+cPCsDh35P+rt7Q3oEGcghv65lHn/g4pLHv8/lzxnz+r4ffeGNdZQCeptFwAAECImfOHULAF/5wMAACAYrHwAAGABZpz5MAvJBwAAVsC2CwAAQHiw8gEAgBUEue0STSsfJB8AAFgB2y4AAADhwcoHAABWEEMrHyQfAABYQCy9asu2CwAAiCiSDwAAEFFsuwAAYAWc+QAAAJHEmQ8AAIAwYeUDAACriKLVi2CQfAAAYAUxdOaDbRcAABBRrHwAAGABHDiNoLq6OmVmZio5OVkOh0OHDx8etX1/f782bdqkjIwMJSUl6dJLL9XOnTsjFC0AAGFihKBECVNXPhobG1VRUaG6ujotXrxY27ZtU2FhoY4dO6Y5c+b47bNixQqdPHlSO3bs0GWXXaauri4NDg5GOHIAADBepiYfNTU1KikpUWlpqSSptrZWL730kurr61VdXT2s/YsvvqiDBw/q/fff19SpUyVJc+fOHXWM/v5+9ff3e3/u6+sL3QQAAAgRtl0iYGBgQK2trSooKPCpLygoUHNzs98+v/zlL5WTk6Mf//jHuuSSSzRv3jzddddd+uSTT0Ycp7q6Wqmpqd5it9tDOg8AAEKCbZfw6+7ultvtVnp6uk99enq6Ojs7/fZ5//339eqrryo5OVnPPvusuru7VV5erj//+c8jnvuoqqpSZWWl9+e+vj4SEAAATGT62y42m83nZ8MwhtUN8Xg8stls2r17t1JTUyWd27r57ne/q61bt2ry5MnD+iQlJSkpKSn0gQMAEEp85yP80tLSFB8fP2yVo6ura9hqyJCZM2fqkksu8SYekpSVlSXDMPThhx+GNV4AAMJp6MxHMCVamLbykZiYKIfDIafTqW9961veeqfTqVtuucVvn8WLF+vpp5/W6dOndcEFF0iS3nnnHcXFxWn27NkBjZ9w5lMlxJv+pjEARET2NP/b2VbRedqi2+HuTyM3FisfkVFZWaknnnhCO3fu1Ntvv60NGzaovb1dZWVlks6d1yguLva2X7lypaZNm6bbb79dx44d06FDh3T33Xfrjjvu8LvlAgAArMfUMx9FRUXq6enRli1b5HK5lJ2draamJmVkZEiSXC6X2tvbve0vuOACOZ1O3XnnncrJydG0adO0YsUKPfDAA2ZNAQCA0IihlQ/TD5yWl5ervLzc7+8aGhqG1S1YsEBOpzPMUQEAEFl85wMAACBMTF/5AAAAYtsFAABEFtsuAAAAYULyAQCAFZh0t0tdXZ0yMzOVnJwsh8Ohw4cPj9j21Vdf1eLFizVt2jRNnjxZCxYs0E9/+tOAx2TbBQAAKzDhzEdjY6MqKipUV1enxYsXa9u2bSosLNSxY8c0Z86cYe2nTJmidevW6aqrrtKUKVP06quvas2aNZoyZYq+973vjXlcVj4AAIhRNTU1KikpUWlpqbKyslRbWyu73a76+nq/7RctWqRbb71VV155pebOnat//ud/1o033jjqaok/JB8AAFiALQRFOnd7++dLf3+/3/EGBgbU2tqqgoICn/qCggI1NzePKea2tjY1Nzfr2muvDWSqJB8AAFhCiM582O12paamekt1dbXf4bq7u+V2u4dd5pqenj7s0tfzzZ49W0lJScrJydHatWtVWloa0FQ58wEAgAWE6lXbjo4OpaSkeOuTkpJG72ez+fxsGMawuvMdPnxYp0+f1uuvv66NGzfqsssu06233jrmWEk+AACYQFJSUnySj5GkpaUpPj5+2CpHV1fXsNWQ82VmZkqS/uZv/kYnT57U5s2bA0o+2HYBAMAKIvyqbWJiohwOx7D70pxOp/Ly8sYetmGMeK5kJKx8AABgFRH+SmllZaVWrVqlnJwc5ebmavv27Wpvb1dZWZkkqaqqSidOnNCuXbskSVu3btWcOXO0YMECSee++/HII4/ozjvvDGhckg8AAGJUUVGRenp6tGXLFrlcLmVnZ6upqUkZGRmSJJfLpfb2dm97j8ejqqoqHT9+XAkJCbr00kv1ox/9SGvWrAloXJIPAAAswKy7XcrLy1VeXu73dw0NDT4/33nnnQGvcvhD8gEAgBXE0K22HDgFAAARxcoHAAAWYNa2ixlIPgAAsAK2XQAAAMKDlQ8AACyAbRcAABBZMbTtQvIBAIAVxFDywZkPAAAQUax8AABgAZz5AAAAkcW2CwAAQHiw8gEAgAXYDEM2Y/zLF8H0jTSSDwAArIBtFwAAgPBg5QMAAAvgbRcAABBZbLsAAACEh+krH3V1dXr44Yflcrl05ZVXqra2Vvn5+V/Y78iRI7r22muVnZ2to0ePBjxu3Jl+xcWPI2AAiEL5F71jdgij2ntmutkh+BXn7o/YWLG07WLqykdjY6MqKiq0adMmtbW1KT8/X4WFhWpvbx+1X29vr4qLi7VkyZIIRQoAQJgZIShRwtTko6amRiUlJSotLVVWVpZqa2tlt9tVX18/ar81a9Zo5cqVys3NjVCkAACE19DKRzAlWpiWfAwMDKi1tVUFBQU+9QUFBWpubh6x35NPPqn33ntP991335jG6e/vV19fn08BAADmMS356O7ultvtVnp6uk99enq6Ojs7/fZ59913tXHjRu3evVsJCWM7rlJdXa3U1FRvsdvtQccOAEDIse0SOTabzednwzCG1UmS2+3WypUrdf/992vevHlj/vtVVVXq7e31lo6OjqBjBgAgHGJhy0Uy8W2XtLQ0xcfHD1vl6OrqGrYaIkmnTp1SS0uL2tratG7dOkmSx+ORYRhKSEjQ/v37dcMNNwzrl5SUpKSkpPBMAgAABMy05CMxMVEOh0NOp1Pf+ta3vPVOp1O33HLLsPYpKSl68803ferq6ur08ssv65lnnlFmZmbYYwYAIGwM41wJpn+UMPU7H5WVlVq1apVycnKUm5ur7du3q729XWVlZZLObZmcOHFCu3btUlxcnLKzs336T58+XcnJycPqAQCINrH0nQ9Tk4+ioiL19PRoy5Ytcrlcys7OVlNTkzIyMiRJLpfrC7/5AQAAoovpXzgtLy9XeXm53981NDSM2nfz5s3avHlz6IMCACDSYuhuF9OTDwAAINk850ow/aOF6a/aAgCA2MLKBwAAVsC2CwAAiCTedgEAAJEVQ9/54MwHAAAxrK6uTpmZmUpOTpbD4dDhw4dHbLtv3z4tW7ZMF198sVJSUpSbm6uXXnop4DFJPgAAsIBg7nUZ75ZNY2OjKioqtGnTJrW1tSk/P1+FhYUjfmPr0KFDWrZsmZqamtTa2qrrr79ey5cvV1tbW0DjknwAAGAFJtxqW1NTo5KSEpWWliorK0u1tbWy2+2qr6/32762tlb33HOPvva1r+nyyy/Xgw8+qMsvv1wvvPBCQOOSfAAAMIH09fX5lP7+fr/tBgYG1NraqoKCAp/6goICNTc3j2ksj8ejU6dOaerUqQHFSPIBAIAFhGrbxW63KzU11Vuqq6v9jtfd3S232z3sJvn09PRhN86P5Cc/+YnOnDmjFStWBDRX3nYBAMAKQvS2S0dHh1JSUrzVSUlJo3az2Wzn/RljWJ0/e/bs0ebNm/X8889r+vTpAYVK8gEAwASSkpLik3yMJC0tTfHx8cNWObq6uoathpyvsbFRJSUlevrpp7V06dKAY2TbBQAAC4j02y6JiYlyOBxyOp0+9U6nU3l5eSP227Nnj2677TY99dRTuvnmm8czVVY+AACwBBM+r15ZWalVq1YpJydHubm52r59u9rb21VWViZJqqqq0okTJ7Rr1y5J5xKP4uJiPfroo7r66qu9qyaTJ09WamrqmMcl+QAAIEYVFRWpp6dHW7ZskcvlUnZ2tpqampSRkSFJcrlcPt/82LZtmwYHB7V27VqtXbvWW7969Wo1NDSMeVySDwAALMCsu13Ky8tVXl7u93fnJxQHDhwY3yDnIfkAAMAKPMa5Ekz/KEHyAQCAFZhw5sMsvO0CAAAiipUPAAAswKYgz3yELJLwI/kAAMAKQvSF02jAtgsAAIgoVj4AALAAs161NQPJBwAAVsDbLgAAAOHBygcAABZgMwzZgjg0GkzfSIvZ5OOTOalKSEg2OwwAE8SpOYlmhzCqpV/6o9khjOq/Mv7e7BD8Ghw8K70bocE8n5Vg+kcJtl0AAEBExezKBwAAVsK2CwAAiKwYetuF5AMAACvgC6cAAADhwcoHAAAWEEtfODV95aOurk6ZmZlKTk6Ww+HQ4cOHR2y7b98+LVu2TBdffLFSUlKUm5url156KYLRAgAQJkPbLsGUKGFq8tHY2KiKigpt2rRJbW1tys/PV2Fhodrb2/22P3TokJYtW6ampia1trbq+uuv1/Lly9XW1hbhyAEAwHiZuu1SU1OjkpISlZaWSpJqa2v10ksvqb6+XtXV1cPa19bW+vz84IMP6vnnn9cLL7ygRYsWRSJkAADCwuY5V4LpHy1MW/kYGBhQa2urCgoKfOoLCgrU3Nw8pr/h8Xh06tQpTZ06dcQ2/f396uvr8ykAAFgO2y7h193dLbfbrfT0dJ/69PR0dXZ2julv/OQnP9GZM2e0YsWKEdtUV1crNTXVW+x2e1BxAwCA4Jh+4NRms/n8bBjGsDp/9uzZo82bN6uxsVHTp08fsV1VVZV6e3u9paOjI+iYAQAIOSMEJUqYduYjLS1N8fHxw1Y5urq6hq2GnK+xsVElJSV6+umntXTp0lHbJiUlKSkpKeh4AQAIp1j6vLppKx+JiYlyOBxyOp0+9U6nU3l5eSP227Nnj2677TY99dRTuvnmm8MdJgAACDFT33aprKzUqlWrlJOTo9zcXG3fvl3t7e0qKyuTdG7L5MSJE9q1a5ekc4lHcXGxHn30UV199dXeVZPJkycrNTXVtHkAABC0GPq8uqnJR1FRkXp6erRlyxa5XC5lZ2erqalJGRkZkiSXy+XzzY9t27ZpcHBQa9eu1dq1a731q1evVkNDQ6TDBwAgdAxJwbwuGz25h/mfVy8vL1d5ebnf352fUBw4cCD8AQEAYALOfAAAAISJ6SsfAABAn70uG8yZj5BFEnYkHwAAWEEMHThl2wUAAEQUKx8AAFiBR9IXf+B79P5RguQDAAAL4G0XAACAMCH5AADACoYOnAZTxqGurk6ZmZlKTk6Ww+HQ4cOHR2zrcrm0cuVKzZ8/X3FxcaqoqBjXmCQfAABYgQnJR2NjoyoqKrRp0ya1tbUpPz9fhYWFPl8X/7z+/n5dfPHF2rRpkxYuXDjuqZJ8AAAwgfT19fmU/v7+EdvW1NSopKREpaWlysrKUm1trex2u+rr6/22nzt3rh599FEVFxcHdacayQcAAFYQopUPu92u1NRUb6murvY73MDAgFpbW1VQUOBTX1BQoObm5rBOlbddAACwghC9atvR0aGUlBRvdVJSkt/m3d3dcrvdSk9P96lPT0/33hofLiQfAABYQKhetU1JSfFJPr6wn8034zEMY1hdqLHtAgBADEpLS1N8fPywVY6urq5hqyGhRvIBAIAVRPhtl8TERDkcDjmdTp96p9OpvLy8UM5sGLZdAACwAo8h2YL4Sqkn8L6VlZVatWqVcnJylJubq+3bt6u9vV1lZWWSpKqqKp04cUK7du3y9jl69Kgk6fTp0/qf//kfHT16VImJibriiivGPC7JBwAAMaqoqEg9PT3asmWLXC6XsrOz1dTUpIyMDEnnPip2/jc/Fi1a5P2/W1tb9dRTTykjI0N/+tOfxjwuyQcAAFYQxFdKvf3Hoby8XOXl5X5/19DQ4GeY4O+QIfkAAMASgkw+FD0Xy8Vs8jHp/2tTgm2S2WFEJVt8vNkhRK24yZPNDmFUtsnJZocwsqkXmR3BqP562XSzQxjVnIQLzQ5hVMnNvzc7BL8GjQGzQ5iQYjb5AADAUkzadjEDyQcAAFbgMRTU1sk43nYxC9/5AAAAEcXKBwAAVmB4zpVg+kcJkg8AAKyAMx8AACCiOPMBAAAQHqx8AABgBWy7AACAiDIUZPIRskjCjm0XAAAQUax8AABgBWy7AACAiPJ4JAXxrQ5P9Hzng20XAAAQUax8AABgBTG07WL6ykddXZ0yMzOVnJwsh8Ohw4cPj9r+4MGDcjgcSk5O1le+8hU9/vjjEYoUAIAwGko+gilRwtTko7GxURUVFdq0aZPa2tqUn5+vwsJCtbe3+21//Phx3XTTTcrPz1dbW5vuvfderV+/Xnv37o1w5AAAYLxMTT5qampUUlKi0tJSZWVlqba2Vna7XfX19X7bP/7445ozZ45qa2uVlZWl0tJS3XHHHXrkkUdGHKO/v199fX0+BQAAy/EYwZcoYVryMTAwoNbWVhUUFPjUFxQUqLm52W+f1157bVj7G2+8US0tLfr000/99qmurlZqaqq32O320EwAAIAQMgxP0CVamJZ8dHd3y+12Kz093ac+PT1dnZ2dfvt0dnb6bT84OKju7m6/faqqqtTb2+stHR0doZkAAAChZAS56hFFZz5Mf9vFZrP5/GwYxrC6L2rvr35IUlKSkpKSgowSAACEimnJR1pamuLj44etcnR1dQ1b3RgyY8YMv+0TEhI0bdq0sMUKAEDYGYaCuqAlilY+TNt2SUxMlMPhkNPp9Kl3Op3Ky8vz2yc3N3dY+/379ysnJ0eTJk0KW6wAAISdxxN8iRKmvu1SWVmpJ554Qjt37tTbb7+tDRs2qL29XWVlZZLOndcoLi72ti8rK9MHH3ygyspKvf3229q5c6d27Nihu+66y6wpAACAAJl65qOoqEg9PT3asmWLXC6XsrOz1dTUpIyMDEmSy+Xy+eZHZmammpqatGHDBm3dulWzZs3Sz372M33nO98xawoAAIRGDG27mH7gtLy8XOXl5X5/19DQMKzu2muv1W9+85swRwUAQGQZHo8M2/i3TnjVFgAAYASmr3wAAACx7QIAACLMY0i22Eg+2HYBAAARxcoHAABWYBiSgjg0GkUrHyQfAABYgOExZASx7WKQfAAAgIAYHgW38sGrtgAAAH6x8gEAgAWw7QIAACIrhrZdYi75GMoMB/VpUN9yiWW2KPo3uNXEGfFmhzAqm8fCO7HufrMjGJV74KzZIYyq75Tb7BBGNWgMmB2CX4PGp5Iis6oQ7D+XBvVp6IIJM5sRTes0IfDhhx/KbrebHQYAIIp0dHRo9uzZYfnbZ8+eVWZmpjo7O4P+WzNmzNDx48eVnJwcgsjCJ+aSD4/Ho48++kgXXnihbDZb0H+vr69PdrtdHR0dSklJCUGE1hZL842luUqxNd9YmqsUW/MN9VwNw9CpU6c0a9YsxcWFb2Xw7NmzGhgIfvUnMTHR8omHFIPbLnFxcWHJXlNSUib8f6g/L5bmG0tzlWJrvrE0Vym25hvKuaampobk74wmOTk5KpKGULHwBi8AAJiISD4AAEBEkXwEKSkpSffdd5+SkpLMDiUiYmm+sTRXKbbmG0tzlWJrvrE012gWcwdOAQCAuVj5AAAAEUXyAQAAIorkAwAARBTJBwAAiCiSjzGoq6tTZmamkpOT5XA4dPjw4VHbHzx4UA6HQ8nJyfrKV76ixx9/PEKRhkYg8z1w4IBsNtuw8vvf/z6CEY/PoUOHtHz5cs2aNUs2m03PPffcF/aJ1mcb6Fyj+blWV1fra1/7mi688EJNnz5d3/zmN/WHP/zhC/tF67Mdz3yj9fnW19frqquu8n5ALDc3V7/+9a9H7ROtz3WiI/n4Ao2NjaqoqNCmTZvU1tam/Px8FRYWqr293W/748eP66abblJ+fr7a2tp07733av369dq7d2+EIx+fQOc75A9/+INcLpe3XH755RGKePzOnDmjhQsX6rHHHhtT+2h+toHOdUg0PteDBw9q7dq1ev311+V0OjU4OKiCggKdOXNmxD7R/GzHM98h0fZ8Z8+erR/96EdqaWlRS0uLbrjhBt1yyy166623/LaP5uc64RkY1de//nWjrKzMp27BggXGxo0b/ba/5557jAULFvjUrVmzxrj66qvDFmMoBTrfV155xZBk/OUvf4lAdOEjyXj22WdHbRPtz3bIWOY6UZ6rYRhGV1eXIck4ePDgiG0myrM1jLHNdyI93y9/+cvGE0884fd3E+m5TjSsfIxiYGBAra2tKigo8KkvKChQc3Oz3z6vvfbasPY33nijWlpa9Omn1r7ueDzzHbJo0SLNnDlTS5Ys0SuvvBLOME0Tzc92vCbCc+3t7ZUkTZ06dcQ2E+nZjmW+Q6L5+brdbv3iF7/QmTNnlJub67fNRHquEw3Jxyi6u7vldruVnp7uU5+enj7i1cednZ1+2w8ODqq7uztssYbCeOY7c+ZMbd++XXv37tW+ffs0f/58LVmyRIcOHYpEyBEVzc82UBPluRqGocrKSl1zzTXKzs4esd1EebZjnW80P98333xTF1xwgZKSklRWVqZnn31WV1xxhd+2E+W5TkQxd6vteNhsNp+fDcMYVvdF7f3VW1Ug850/f77mz5/v/Tk3N1cdHR165JFH9I1vfCOscZoh2p/tWE2U57pu3Tr97ne/06uvvvqFbSfCsx3rfKP5+c6fP19Hjx7VX//6V+3du1erV6/WwYMHR0xAJsJznYhY+RhFWlqa4uPjh/23/q6urmHZ9JAZM2b4bZ+QkKBp06aFLdZQGM98/bn66qv17rvvhjo800Xzsw2FaHuud955p375y1/qlVde0ezZs0dtOxGebSDz9Sdanm9iYqIuu+wy5eTkqLq6WgsXLtSjjz7qt+1EeK4TFcnHKBITE+VwOOR0On3qnU6n8vLy/PbJzc0d1n7//v3KycnRpEmTwhZrKIxnvv60tbVp5syZoQ7PdNH8bEMhWp6rYRhat26d9u3bp5dfflmZmZlf2Cean+145utPtDzf8xmGof7+fr+/i+bnOuGZdNA1avziF78wJk2aZOzYscM4duyYUVFRYUyZMsX405/+ZBiGYWzcuNFYtWqVt/37779vfOlLXzI2bNhgHDt2zNixY4cxadIk45lnnjFrCgEJdL4//elPjWeffdZ45513jP/+7/82Nm7caEgy9u7da9YUxuzUqVNGW1ub0dbWZkgyampqjLa2NuODDz4wDGNiPdtA5xrNz/Vf//VfjdTUVOPAgQOGy+Xylo8//tjbZiI92/HMN1qfb1VVlXHo0CHj+PHjxu9+9zvj3nvvNeLi4oz9+/cbhjGxnutER/IxBlu3bjUyMjKMxMRE46tf/arPK2yrV682rr32Wp/2Bw4cMBYtWmQkJiYac+fONerr6yMccXACme9DDz1kXHrppUZycrLx5S9/2bjmmmuMX/3qVyZEHbih1w3PL6tXrzYMY2I920DnGs3P1d88JRlPPvmkt81EerbjmW+0Pt877rjD+/+bLr74YmPJkiXexMMwJtZznehshvHZ6RsAAIAI4MwHAACIKJIPAAAQUSQfAAAgokg+AABARJF8AACAiCL5AAAAEUXyAQAAIorkAwAARBTJBwAAiCiSDyCGXHfddaqoqBhW/9xzz3HFOICIIfkAAAARRfIBwMdvf/tbXX/99brwwguVkpIih8OhlpYWs8MCMIEkmB0AAGv5p3/6Jy1atEj19fWKj4/X0aNHNWnSJLPDAjCBkHwA8NHe3q67775bCxYskCRdfvnlJkcEYKJh2wWAj8rKSpWWlmrp0qX60Y9+pPfee8/skABMMCQfQAxJSUlRb2/vsPq//vWvSklJkSRt3rxZb731lm6++Wa9/PLLuuKKK/Tss89GOlQAExjJBxBDFixY4Pfw6BtvvKH58+d7f543b542bNig/fv369vf/raefPLJSIYJYIIj+QBiSHl5ud577z2tXbtWv/3tb/XOO+9o69at2rFjh+6++2598sknWrdunQ4cOKAPPvhAR44c0RtvvKGsrCyzQwcwgdgMwzDMDgJA5LS2tmrTpk1qa2vT2bNnNW/ePH3/+9/XP/7jP2pgYECrV6/WkSNHdPLkSaWlpenb3/62Hn74YSUnJ5sdOoAJguQDAABEFNsuAAAgokg+AABARJF8AACAiCL5AAAAEUXyAQAAIorkAwAARBTJBwAAiCiSDwAAEFEkHwAAIKJIPgAAQESRfAAAgIj6vx+IJLGv3NdkAAAAAElFTkSuQmCC",
+      "text/plain": [
+       "<Figure size 640x480 with 2 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "gap_da.plot(x=\"Us\", y=\"Vs\")"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 8,
@@ -225,7 +411,7 @@
     }
    ],
    "source": [
-    "gap_da.plot(x='Us', y='Vs')"
+    "gap_da.plot(x=\"Us\", y=\"Vs\")"
    ]
   },
   {
@@ -235,7 +421,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "gap_da.to_netcdf('./data/graphene_example.nc')"
+    "gap_da.to_netcdf(\"./data/graphene_example.nc\")"
    ]
   }
  ],
@@ -255,7 +441,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.11.6"
   },
   "widgets": {
    "application/vnd.jupyter.widget-state+json": {
-- 
GitLab