RuCl documentation update

notebooks/RuClExample.ipynb

@@ -2,11 +2,37 @@
  "cells": [
   {
    "cell_type": "markdown",
-   "id": "c375a2e9-f691-4aba-8f0f-ac466ce86f33",
+   "id": "a6665d70-aa98-4284-a53c-e315a466ed95",
    "metadata": {},
    "source": [
     "# RuCl Full Scale Example\n",
-    "We can now implement a full scale example, with details described here https://www.nature.com/articles/s41535-019-0203-y"
+    "The search for Kitaev Spin Liquid Materials is a very active area of Material Science Research.  One example of a candidate material that demonstrates promise as a Kitaev Spin Liquid is $\\alpha$-RuCl$_3$. In this notebook, we explore the necessary steps to simulate the dynamics of this material on a Fault Tolerant Quantum Computer, assuming a gate-based architecture. For a good introduction to the properties of the material see the review article [here](https://www.nature.com/articles/s41535-019-0203-y).  The Hamiltonian for $\\alpha$-RuCl$_3$ can be most generically represented in the form\n",
+    "\\begin{equation}\n",
+    "\\begin{split}\n",
+    "  H_{material} &= K_{x} \\sum_{ij} S^{x}_{i} S^{x}_{j} + K_{y} \\sum_{ij} S^{y}_{i} S^{y}_{j} + K_{z} \\sum_{ij} S^{z}_{i} S^{z}_{j} + J \\sum_{ij} \\bf{S_i} \\cdot \\bf{S_j}\\\\\n",
+    "    &+ \\Gamma_z \\sum_{ij} (S^{x}_{i} S^{y}_{j} + S^{y}_{i} S^{x}_{j}) + \\Gamma_y \\sum_{ij} (S^{z}_{i} S^{x}_{j} + S^{x}_{i} S^{z}_{j}) + \\Gamma_x \\sum_{ij} (S^{y}_{i} S^{z}_{j} + S^{z}_{i} S^{y}_{j}) \\\\\n",
+    "    &+ \\Gamma_z' \\sum_{ij} (S^{x}_{i} S^{z}_{j} + S^{z}_{i} S^{x}_{j}) + (S^{x}_{i} S^{y}_{j} + S^{y}_{i} S^{x}_{j}) \\\\\n",
+    "    &+ \\Gamma_y' \\sum_{ij} (S^{y}_{i} S^{x}_{j} + S^{x}_{i} S^{y}_{j}) + (S^{y}_{i} S^{z}_{j} + S^{z}_{i} S^{y}_{j}) \\\\\n",
+    "    &+ \\Gamma_x' \\sum_{ij} (S^{z}_{i} S^{x}_{j} + S^{x}_{i} S^{z}_{j}) + (S^{z}_{i} S^{y}_{j} + S^{y}_{i} S^{z}_{j}) \\\\\n",
+    "    &+ A \\sum_{i} (S^{z}_i)^{2} + \\sum_{i} \\bf{S_i \\cdot H_i}\n",
+    "\\end{split}\n",
+    "\\end{equation}\n",
+    "\n",
+    "The terms $S^{x}_{i}$, $S^{y}_{i}$, and $S^{z}_{i}$ represent the Pauli operators acting on site $i$.  The terms $K_{x}$, $K_{y}$, and $K_{z}$ represent the strength of the Kitaev interaction between two sites in a given direction.  The bold terms $\\bf{S_i}$ $ = [S^x_i, S^y_i, S^z_i]$ are useful for defining the Heisenberg interaction terms with strength $J$.  The terms $\\Gamma$ represent the strength of off-diagonal symmetric exchange interactions between nearest neighboring sites, and the terms $\\Gamma'$ represent the effect of trigonal distortion.  Lastly, the term $A$ represents the effect of single-ion anisotropy.  Finally, $\\bf{H}$ is a vector representing the strength of the Zeeman energy.\n",
+    "\n",
+    "The connectivity of the Hamiltonian is defined directionally, as shown in Figure 1 (obtained from [[1]](https://doi.org/10.1038/s41535-019-0203-y)).  As a result of the experiment being performed on the material, we need to include a time-varying Hamiltonian component, corresponding to the time-varying transverse field of the magnet operating on the material. This Hamiltonian is\n",
+    "\\begin{equation}\n",
+    "  H_{field}(t) = t \\sum_i S^x_i\n",
+    "\\end{equation}\n",
+    "\n",
+    "We can then construct the complete time-varying Hamiltonian we want to simulate with \n",
+    "\\begin{equation}\n",
+    "  H(t) = H_{material} + H_{field}(t)\n",
+    "\\end{equation}\n",
+    "\n",
+    "##### Figure 1\n",
+    "\n",
+    "![RuCl_Lattice](EmbeddedFigures/RuclLattice.jpg)"
    ]
   },
   {
@@ -16,16 +42,7 @@
    "metadata": {
     "tags": []
    },
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/Users/zain/anaconda3/lib/python3.11/site-packages/attr/_make.py:918: RuntimeWarning: Running interpreter doesn't sufficiently support code object introspection.  Some features like bare super() or accessing __class__ will not work with slotted classes.\n",
-      "  set_closure_cell(cell, cls)\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "import os\n",
     "import time\n",
@@ -53,6 +70,32 @@
     "from qca.utils.hamiltonian_utils import flatten_nx_graph, assign_hexagon_labels, pyliqtr_hamiltonian_to_openfermion_qubit_operator"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "38dcb3a9-bdd1-4829-a07e-2528ddce2206",
+   "metadata": {},
+   "source": [
+    "First we will define the functions that we are going to use to estimate the resources for simulating the dynamics of the Hamiltonian, $H$ described above. Since most of the common methods Hamiltonian simulation currently in use, for example Quantum Signal Processing (QSP) and Product formulas (e.g. Trotterization), we need to first break up the Hamiltonian into time independent steps using a Magnus expansion.  Here, we only use the first order Magnus expansion, which ultimately bottlenecks our precision, but we allow the user to specify a number of steps to ensure that the simulation is accurate enough for practical purposes.\n",
+    "\n",
+    "### Magnus Expansion\n",
+    "The operation of the first order Magnus expansion is defined as follows:\n",
+    "\\begin{equation}\n",
+    "U(s, s_0) \\approx U_{mag_1} = \\exp{\\left(-i \\int_{t0}^{t1} H(s) ds \\right)}\n",
+    "\\end{equation}\n",
+    "\n",
+    "implying that we can treat the Hamiltonian as the argument of the exponential:\n",
+    "\\begin{equation}\n",
+    "H(s, s_0) \\approx H_{mag_1} = \\int_{t0}^{t1} H(s) ds\n",
+    "\\end{equation}\n",
+    "\n",
+    "Note that when H has non-commuting component terms, $H = H_1 + H_2$ where $[H_1, H_2] \\neq 0$, this  serves only as a first order method and higher order terms can be introduced for more precision.  For the sake of simplicity we choose to use the First order approximation.  As a further approximating step, in order to reduce overheads of compiling multiple circuits, we currently treat the time varying portion as a constant, representing a value somewhere in the middle of the dynamics.  In the future, we will likely improve these approximations to fully explore the time varying component and the potential improvement of using higher-order magnus expansions.\n",
+    "\n",
+    "### Simulation Algorithms\n",
+    "In this notebook, we use two different algorithms for simulating the dynamics of $H$: QSP and second order Trotterization.  We explore both of these methods to better compare the resource requirements for each algorithm and see where tradeoffs, such as the parallelizability of T gates, exist.  It should be noted that the bounds which we use for estimating the Trotter error are loose, and it has been shown that using a lower number of steps is often sufficient.  For the implementation of QSP, we use [pyLIQTR](https://github.com/isi-usc-edu/pyLIQTR) and for the implementation of Trotterization, we use [openfermion](https://quantumai.google/openfermion).  Because the term ordering used by default is almost antagonistic with respect to circuit depth for most Hamiltonians with a lattice structure, we also define our own term ordering based on edge coloring.  While this may not be perfectly optimal, this coloring approach tends to get rather close empirically.\n",
+    "\n",
+    "Due to the large circuit sizes required for implementing both of these algorithms, we take approaches to reduce the portion of the circuit directly stored in memory at any given time, and extrapolate the results for various subcircuits to obtain an approximate final estimate.  In the case of QSP, the circuits are provided in the form of high level circuit abstractions which may be decomposed later.  There is a high degree of repetition of the blocks in this high level abstraction, so we take one instance of each block, decompose it to its required clifford + T operations, and multiply that count by the number of occurences of the block.  For Trotterization, since each step is iterated multiple times, we simply compute the resources required for a single step and multiply those resource estimates by the number of steps."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 2,
@@ -63,8 +106,8 @@
    "outputs": [],
    "source": [
     "t_init = time.perf_counter()\n",
-    "def estimate_qsp(pyliqtr_hamiltonian, timesteps, energy_precision, outdir, hamiltonian_name=\"hamiltonian\", write_circuits=False):\n",
-    "    timestep_of_interest=1 #for magnus like argument\n",
+    "def estimate_qsp(pyliqtr_hamiltonian, timesteps, energy_precision, outdir, hamiltonian_name=\"hamiltonian\", write_circuits=False, numsteps=1000):\n",
+    "    timestep_of_interest=timesteps/numsteps #for magnus like argument\n",
     "    t0 = time.perf_counter()\n",
     "    random.seed()\n",
     "    angles_response = Angler_fourier_response(tau=timestep_of_interest*pyliqtr_hamiltonian.alpha, eps=energy_precision, random=True, silent=True)\n",
@@ -262,6 +305,14 @@
     "    return cpt_trotter"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "3263215d-b585-4728-a5f3-96baf4faf050",
+   "metadata": {},
+   "source": [
+    "Next, we need to define the directional labels for the hexagonal lattice and add next-nearest and next-next-nearest neighbor interactions to the graph so we can construct the Hamiltonian from the connectivity graph later on.  This should match the connectivity defined in figure 1 above.  Note that we only consider cross-plaquette next-next-nearest neighbor interactions.  For the sake of visualization, we use a 3x3 lattice.  For the actual simulations, we will use a much larger lattice, defined below."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 3,
@@ -314,6 +365,14 @@
     "nx.draw_networkx_edge_labels(g_rucl, pos,edge_labels = edge_labels);"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "3e1528a0-de81-4807-820f-d273871e7cfe",
+   "metadata": {},
+   "source": [
+    "With the connectivity defined, we can construct a variety of Hamiltonians using the terms defined in [[1]](https://doi.org/10.1038/s41535-019-0203-y).  Some of the models contain more terms and thus take longer to run.  For the purposes of demonstration, row 13 tends to have the shortest runtime, so we use that for demonstration purposes here.  To run all rows, simply loop over all indices rather than just the 13th index.  The number after the terms in the following graph indicate how far away the interactions are.  For example J1 represents the nearest neighbor Heisenberg interactions, J2 represents the next-nearest neighbor Heisenberg interactions, and J3 represents the next-next-nearest neighbor Heisenberg interactions.  Any terms not present in the table that appear in the generic Hamiltonian defined above are treated as 0."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 4,
@@ -683,6 +742,17 @@
     "df_rucl"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "c8f84ab3-4a8d-4893-886c-89bcfb41a1b8",
+   "metadata": {},
+   "source": [
+    "##### Table 1\n",
+    "\n",
+    "\n",
+    "Now we may construct the Hamiltonian in a form usable by pyLIQTR using the values defined in the table above.  We also allow for an additional time varying portion to represent the longitudinal field term.  In the future, we would like to extend the capability of this notebook to further explore implementing this term in a more precise manner with higher order Magnus expansions and selecting multiple time-slices for a more accurate approximation."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 5,
@@ -811,6 +881,14 @@
     "    return H"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "7969824d-de0c-4e29-9572-89d6fa55612d",
+   "metadata": {},
+   "source": [
+    "Finally, it is time to implement the circuit for simulating the time dynamics of $H$, parameterized by the rows of Table 1.  First we perform the simulation using the second order Trotter expansion.  We choose an energy precision of 1e-3, which is believed to be sufficiently low to allow for lower circuit depths while still getting empirically accurate results.  We construct a lattice of size 32x32 honeycomb cells to hopefully have a sufficiently large mass to mitigate finite size effects.  This is believed to be the lower end of what lattice size would be useful, and larger may be better in some cases.  Finally, we plot the a histogram of the T-widths that occur on each layer of the circuit.  The x-axis representing the number of parallel T-gates and the y-axis representing the number of layers where that T-width is found.  This is done to determine the value of the parallelizability of T factories on a theoretical Fault Tolerant Device using a Clifford + T gateset."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 6,
@@ -859,7 +937,7 @@
     "        os.makedirs(widthdir)\n",
     "\n",
     "    timesteps=1000\n",
-    "    H_rucl = generate_rucl_hamiltonian(32, df_rucl.iloc[i])\n",
+    "    H_rucl = generate_rucl_hamiltonian(32, df_rucl.iloc[i], field_x=lambda s: 1)\n",
     "    H_rucl_pyliqtr = pyH(H_rucl)\n",
     "    openfermion_hamiltonian_rucl = pyliqtr_hamiltonian_to_openfermion_qubit_operator(H_rucl_pyliqtr)\n",
     "    \n",
@@ -877,6 +955,14 @@
     "                   0)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "fee1aaac-e081-4008-bd46-45f3051aebd6",
+   "metadata": {},
+   "source": [
+    "Finally we perform the same experiment for the Quantum Signal Processing algorithm.  We omit plotting in this case because we perform estimation of the circuit by decomposing sub-circuits of the overall algorithm.  This results in some sub-circuits where parallelizability may appear deceptively high, or no T gates may even be present.  Emperically, for the portions of the circuit which compose the bulk of the circuit (SELECT and REFLECT operators), the T width tends to be near 1."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 7,
@@ -922,7 +1008,7 @@
     "    if not os.path.exists(widthdir):\n",
     "        os.makedirs(widthdir)\n",
     "    timesteps=1000\n",
-    "    H_rucl = generate_rucl_hamiltonian(32, df_rucl.iloc[i])\n",
+    "    H_rucl = generate_rucl_hamiltonian(32, df_rucl.iloc[i], field_x=lambda s: 1)\n",
     "    H_rucl_pyliqtr = pyH(H_rucl)\n",
     "    \n",
     "    print(\"Estimating RuCl row \" + str(i) + \" using QSP\")\n",
@@ -956,18 +1042,27 @@
     "print(\"Total time to run notebook (only using the fastest row): \" + str(t_end - t_init))"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "d231d98f-920e-42df-8b7d-27b4f2ff46b0",
+   "metadata": {},
+   "source": [
+    "### References\n",
+    "[1] Laurell, P., Okamoto, S. Dynamical and thermal magnetic properties of the Kitaev spin liquid candidate α-RuCl3. npj Quantum Mater. 5, 2 (2020). https://doi.org/10.1038/s41535-019-0203-y"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "09c39464-1b86-4369-8f06-b12f4db46772",
+   "id": "4b08c06a-c632-41f9-a4cd-843b0108599c",
    "metadata": {},
    "outputs": [],
    "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "base",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -981,7 +1076,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.8.8"
   }
  },
  "nbformat": 4,

pyproject.toml

@@ -13,7 +13,7 @@ dependencies = [
     "openfermion",
     "openfermionpyscf",
     "pandas",
-    "pyLIQTR >=1.0.0"
+    "pyLIQTR ==1.0.0"
 ]
 requires-python = ">=3.9"
 description = "Documentation of Applications for Quantum Computers"