DMRG.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "18ed9f64-62cd-4fe6-b92a-eab78e796661",
   "metadata": {},
   "source": [
    "See https://kaihsinwu.gitlab.io/Cytnx_doc/index.html for more information of the Cytnx library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b1b891ef-1e3b-494e-be82-4529dae8ba2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"1\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a88f1cb0-eb8a-45fd-a395-1e4aa2e9147b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MKL interface code: 1 >> using [ilp64] interface\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import cytnx\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb6a0add-cbed-4b6e-bfde-74ef7b137727",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Density matrix renormalization group"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47430e43-ae21-49df-88ec-677cc85bf882",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1. Define the matrix product operator (MPO) for the Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9667d28-e00c-4137-8564-10f3532f96b9",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1.1 Define operators"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76f00a32-2181-4b2f-a0d8-3a25b64bf64e",
   "metadata": {},
   "source": [
    "#### $$ \\hat{S}^{+}=\\left[\\begin{array}{cc}\n",
    "0 & 1\\\\\n",
    "0 & 0\n",
    "\\end{array}\\right]\n",
    ", \\qquad\n",
    "\\hat{S}^{-}=\\left[\\begin{array}{cc}\n",
    "0 & 0\\\\\n",
    "1 & 0\n",
    "\\end{array}\\right]$$\n",
    "#### $$ \\hat{S}^{z}=\\frac{1}{2}\\left[\\begin{array}{cc}\n",
    "1 & 0\\\\\n",
    "0 & -1\n",
    "\\end{array}\\right]\n",
    ", \\qquad\n",
    "\\hat{S}^{x}=\\frac{1}{2}\\left[\\begin{array}{cc}\n",
    "0 & 1\\\\\n",
    "1 & 0\n",
    "\\end{array}\\right]\n",
    ", \\qquad\n",
    "\\hat{I}=\\left[\\begin{array}{cc}\n",
    "1 & 0\\\\\n",
    "0 & 1\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb689c24-a1d9-40e6-adb8-e7d723dd488e",
   "metadata": {},
   "source": [
    "<!-- ![alt text](szsx.png \"Title\") -->\n",
    "<img src=\"figures/szsx.png\" alt=\"drawing\" width=\"300\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d33cb8be-1126-4d1e-bd1b-1e6c97cfeff7",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 2  # Physical dimension\n",
    "\n",
    "# The operators are define in type <Tensor> (rather than <UniTensor>)\n",
    "Sx = cytnx.physics.spin(0.5,'x').real()      # Sx\n",
    "Sz = cytnx.physics.spin(0.5,'z').real()      # Sz\n",
    "I = cytnx.eye(2)                               # Identity operator\n",
    "\n",
    "# S+\n",
    "Sp = cytnx.zeros([2,2])\n",
    "Sp[0,1] = 1\n",
    "\n",
    "# S-\n",
    "Sm = cytnx.zeros([2,2])\n",
    "Sm[1,0] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e23a764-7e2b-4e37-83c1-78cc685ac5d4",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1.2 Define Hamiltonian MPO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bd61770-3366-4354-a480-3d990d3a05fe",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Transverse-field Ising model\n",
    "#### $$\\hat{H} = -\\sum_{i=1}^{N-1} \\hat{S}^x_i\\hat{S}^x_{i+1} - h_z\\sum_{i=1}^N \\hat{S}^z_i$$\n",
    "#### MPO tensor\n",
    "#### $$\\hat{H} = L_0 M_1 M_2 \\cdots M_N R_0,$$\n",
    "#### $$\n",
    "L_{0}=\\left[\\begin{array}{ccc}\n",
    "0 & 0 & 1\\end{array}\\right]\n",
    ", \\qquad\n",
    "M_i=\\left[\\begin{array}{ccc}\n",
    "\\hat{I_i}\\\\\n",
    "\\hat{S}^{x}_i & 0\\\\\n",
    "-h_{z}\\hat{S}^{z}_i & -\\hat{S}^{x}_i & \\hat{I}_i\n",
    "\\end{array}\\right]\n",
    ", \\qquad\n",
    "R_{0}=\\left[\\begin{array}{c}\n",
    "1\\\\\n",
    "0\\\\\n",
    "0\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c93a5cba-bfd7-4391-8155-b0e6b8c4b865",
   "metadata": {},
   "source": [
    "<img src=\"figures/mpo.png\" alt=\"drawing\" width=\"600\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "55ef9566-fb3e-4775-99cf-f9a3d3803a81",
   "metadata": {},
   "outputs": [],
   "source": [
    "dM = 3               # MPO bond dimension\n",
    "def TFIM_MPO_tensor (hz, hx):\n",
    "    M = cytnx.zeros([dM, dM, d, d])\n",
    "    M[0,0] = I\n",
    "    M[2,2] = I\n",
    "    M[1,0] = Sx\n",
    "    M[2,1] = -Sx\n",
    "    M[2,0] = -hz*Sz - hx*Sx\n",
    "    M = cytnx.UniTensor (M, rowrank=2)\n",
    "    return M\n",
    "\n",
    "L0 = cytnx.UniTensor (cytnx.zeros([dM,1,1]), rowrank=1)  # Left boundary\n",
    "R0 = cytnx.UniTensor (cytnx.zeros([dM,1,1]), rowrank=1)  # Right boundary\n",
    "L0.get_block_()[dM-1,0,0] = 1.\n",
    "R0.get_block_()[0,0,0] = 1."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90c5081c-2874-4f46-8c6f-1524bf5ce3d7",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Initialize MPS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dac0266-7ba3-4f31-a252-f5bcd8b760c2",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 2.1 Initialize MPS as a product state\n",
    "### $$ |\\psi\\rangle = |\\uparrow\\downarrow\\uparrow\\downarrow\\cdots\\uparrow\\downarrow\\rangle$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41ad93f0-6763-4b49-8c0e-0d2069fcbe51",
   "metadata": {},
   "source": [
    "<img src=\"figures/mps.png\" alt=\"drawing\" width=\"450\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a058de41-00e6-4ef1-aa58-b5d81a9e4632",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize MPS as a product state\n",
    "def get_init_MPS (Nsites):\n",
    "    psi = [None for i in range(Nsites)]\n",
    "    for i in range(Nsites):\n",
    "        A = cytnx.zeros([1,d,1])\n",
    "        spin = (i % 2)  # 0 for spin up and 1 for spin down\n",
    "        A[0,spin,0] = 1\n",
    "        #A[0,0,0] = 0.5**0.5\n",
    "        #A[0,1,0] = 0.5**0.5\n",
    "        psi[i] = cytnx.UniTensor (A, rowrank=2)\n",
    "        psi[i].set_labels([str(2*i),str(2*i+1),str(2*i+2)])\n",
    "    return psi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ae70ae1-6732-456d-bba4-940dfb6dc917",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 3. Define the effective Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4a00466-0632-4cb8-89cd-c28b9c75d5e5",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.1 Compute the environment tensors LR, and make the MPS **left orthogonal**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7de71c4-8927-4d78-a84e-924bd5bdbce1",
   "metadata": {},
   "source": [
    "<img src=\"figures/leftotho.png\" alt=\"drawing\" width=\"700\"/>\n",
    "<img src=\"figures/envcontract.png\" alt=\"drawing\" width=\"400\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "58f04054-bae8-475a-b6a2-fbbbdfd0928f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#LR[0]:        LR[1]:            LR[2]:\n",
    "#\n",
    "#   -----      -----A[0]---     -----A[1]---\n",
    "#   |          |     |          |     |\n",
    "#  ML----     LR[0]--M-----    LR[1]--M-----      ......\n",
    "#   |          |     |          |     |\n",
    "#   -----      -----A*[0]--     -----A*[1]--\n",
    "#\n",
    "#\n",
    "def get_LR_envir_tensors (psi, M):\n",
    "    Nsites = len(psi)\n",
    "\n",
    "    LR = [None for i in range(Nsites+1)]\n",
    "    LR[0]  = L0\n",
    "    LR[-1] = R0\n",
    "\n",
    "    anet = cytnx.Network()\n",
    "    anet.FromString([\"L: -2,-1,-3\",\n",
    "                     \"A: -1,-4,1\",\n",
    "                     \"M: -2,0,-4,-5\",\n",
    "                     \"Aconj: -3,-5,2\",\n",
    "                     \"TOUT: 0;1,2\"])\n",
    "    for p in range(Nsites - 1):\n",
    "        \n",
    "        # Make A[p] left orthogonal\n",
    "        s, u, vt = cytnx.linalg.Svd(psi[p])\n",
    "\n",
    "        # Setting psi[p] = u\n",
    "        psi[p].put_block(u.get_block())\n",
    "        \n",
    "        # Put the orthogonality center to the right site\n",
    "        r = cytnx.Contract (s, vt)\n",
    "        psi[p+1].put_block(cytnx.Contract (r, psi[p+1]).get_block())\n",
    "        \n",
    "        # Compute the left boundary tensor on site p\n",
    "        anet.PutUniTensors ([\"L\",\"A\",\"Aconj\",\"M\"], [LR[p],psi[p],psi[p].Conj(),M])\n",
    "        LR[p+1] = anet.Launch()\n",
    "    return LR"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9614df74-9989-4f18-aabc-a3ee430819cf",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.2 Define effective Hamiltonian and $\\hat{H}|\\psi\\rangle$ operation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62b32596-5ef6-4788-943d-9e67828af403",
   "metadata": {},
   "source": [
    "<img src=\"figures/heff.png\" alt=\"drawing\" width=\"600\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "72b01134-de05-498a-815e-046c1932710b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# An effective Hamiltonian must:\n",
    "# 1. Inherit <cytnx.LinOp> class\n",
    "# 2. Has a function <matvec> that implements H|psi> operation\n",
    "\n",
    "class eff_Hamilt (cytnx.LinOp):\n",
    "    def __init__ (self, L, M1, M2, R):\n",
    "        # Initialization\n",
    "        cytnx.LinOp.__init__(self,\"mv\", 0)\n",
    "\n",
    "        # Define network for H|psi> operation\n",
    "        anet = cytnx.Network()\n",
    "        anet.FromString([\"psi: -1,-2,-3,-4\",\n",
    "                         \"L: -5,-1,0\",\n",
    "                         \"R: -7,-4,3\",\n",
    "                         \"M1: -5,-6,-2,1\",\n",
    "                         \"M2: -6,-7,-3,2\",\n",
    "                         \"TOUT: 0,1;2,3\"])\n",
    "        anet.PutUniTensors([\"L\",\"M1\",\"M2\",\"R\"],[L,M1,M2,R])\n",
    "        self.anet = anet\n",
    "\n",
    "    # Define H|psi> operation\n",
    "    def matvec (self, v):\n",
    "        lbl = v.labels(); \n",
    "        self.anet.PutUniTensor(\"psi\",v);\n",
    "        out = self.anet.Launch()\n",
    "        out.set_labels(lbl);   # Make sure the input label match output label\n",
    "        return out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e4bd922-1bdd-4082-a001-bceeafe5e3e0",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 4. Find the ground state by DMRG"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b70194ee-2ced-4d84-9a3b-1f1cccee2e51",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4.1 DMRG algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a07c699c-8583-41fb-88e9-9dfacceb98da",
   "metadata": {},
   "source": [
    "<img src=\"figures/dmrg_sweep.png\" alt=\"drawing\" width=\"600\"/>\n",
    "<img src=\"figures/twositecontract.png\" alt=\"drawing\" width=\"650\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "55d773e3-a3fd-42c1-a2d1-6628812e7cf9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy.sparse.linalg import LinearOperator, eigsh, eigs\n",
    "\n",
    "# Optimize phi using Scipy function (Not used)\n",
    "def optimize_psi_scipy(psivec, L, M1, M2, R):\n",
    "    \n",
    "    lbl=psivec.labels()\n",
    "    pshape = [L.shape()[1],M1.shape()[2],M2.shape()[2],R.shape()[1]]\n",
    "    \n",
    "    effH = eff_Hamilt (L, M1, M2, R)\n",
    "    def linop (vec):\n",
    "        vec = vec.reshape(pshape)\n",
    "        vec = cytnx.UniTensor(cytnx.from_numpy(vec), rowrank = 2)\n",
    "        vec = effH.matvec (vec)\n",
    "        return vec.get_block().numpy().flatten()\n",
    "\n",
    "    dim = L.shape()[1] * M1.shape()[2] * M2.shape()[2] * R.shape()[1]\n",
    "    TensorOp = LinearOperator ((dim, dim), matvec=linop, dtype=np.float64)\n",
    "    energy, psivec = eigsh (TensorOp, k=1, which='SA', v0=psivec.get_block().numpy().flatten(),\n",
    "                            ncv=None, maxiter=2, tol=1e-2)\n",
    "    psivec = cytnx.UniTensor(cytnx.from_numpy(psivec.reshape(pshape)))\n",
    "    psivec.set_labels(lbl)\n",
    "    \n",
    "    return psivec, energy\n",
    "\n",
    "def dmrg (psi, M, maxdims, cutoff):\n",
    "    # Define the network for computing the environment tensors\n",
    "    # Network for computing right environment tensors\n",
    "    R_env_net = cytnx.Network()\n",
    "    R_env_net.FromString([\"R: -2,-1,-3\",\n",
    "                          \"B: 1,-4,-1\",\n",
    "                          \"M: 0,-2,-4,-5\",\n",
    "                          \"B_Conj: 2,-5,-3\",\n",
    "                          \"TOUT: 0;1,2\"])\n",
    "    # Network for computing left environment tensors\n",
    "    L_env_net = cytnx.Network()\n",
    "    L_env_net.FromString([\"L: -2,-1,-3\",\\\n",
    "                          \"A: -1,-4,1\",\\\n",
    "                          \"M: -2,0,-4,-5\",\\\n",
    "                          \"A_Conj: -3,-5,2\",\\\n",
    "                          \"TOUT: 0;1,2\"])\n",
    "    anets = [R_env_net, L_env_net]\n",
    "\n",
    "\n",
    "\n",
    "    # Define the links to update for a sweep\n",
    "    # First do a right-to-left and then a left-to-right sweep\n",
    "    ranges = [range(Nsites-2,-1,-1), range(Nsites-1)]\n",
    "\n",
    "    # For printing information\n",
    "    verbose = [\"[r->l]\", \"[l->r]\"]\n",
    "\n",
    "\n",
    "    # Get the environment tensors\n",
    "    LR = get_LR_envir_tensors (psi, M)\n",
    "    \n",
    "    ens = []\n",
    "    for k in range(len(maxdims)):                                                            # For each sweep\n",
    "        chi = maxdims[k]                                                                     # Read bond dimension\n",
    "        print ('Sweep',k,', chi='+str(chi))\n",
    "        for lr in [0,1]:\n",
    "            for p in ranges[lr]:                                                                # For each bond\n",
    "                # Compute the current psi\n",
    "                phi = cytnx.Contract (psi[p], psi[p+1])\n",
    "                phi.set_rowrank(2);\n",
    "                \n",
    "                # Define the effective Hamiltonian\n",
    "                effH = eff_Hamilt (LR[p], M, M, LR[p+2])\n",
    "\n",
    "                # Find the ground state for the current bond\n",
    "                enT, phi = cytnx.linalg.Lanczos(effH, phi, method = \"Gnd\", CvgCrit=999, Maxiter=2)\n",
    "                en = enT.item()                                                                 # Tensor to number\n",
    "                \n",
    "                # phi, en = optimize_psi_scipy (phi, LR[p], M, M, LR[p+2])\n",
    "                \n",
    "                # Store the energy\n",
    "                ens.append(en);\n",
    "\n",
    "                # SVD and truncate the wavefunction psi\n",
    "                s, u, vt = cytnx.linalg.Svd_truncate(phi, keepdim=chi, err=cutoff)      \n",
    "\n",
    "                # Save labels of psi[p] and psi[p+1] for later use\n",
    "                lbl1 = psi[p].labels()\n",
    "                lbl2 = psi[p+1].labels()\n",
    "                \n",
    "                # Setting psi[p] = u, psi[p+1] = vt\n",
    "                psi[p] = u\n",
    "                psi[p+1] = vt\n",
    "                \n",
    "                # Normalize the singular values\n",
    "                s = s/s.get_block_().Norm().item()\n",
    "\n",
    "                if lr == 0:\n",
    "                    # Absorb s into next neighbor\n",
    "                    psi[p] = cytnx.Contract(psi[p],s)\n",
    "\n",
    "                    # Set the network for the right environment tensor on the current site\n",
    "                    anets[lr].PutUniTensors([\"R\",\"B\",\"M\",\"B_Conj\"],[LR[p+2],psi[p+1],M,psi[p+1].Conj()])                \n",
    "                if lr == 1:\n",
    "                    # Absorb s into next neighbor\n",
    "                    psi[p+1] = cytnx.Contract(s,psi[p+1])\n",
    "                    # Set the network for the left environment tensor on the current site\n",
    "                    anets[lr].PutUniTensors([\"L\",\"A\",\"A_Conj\",\"M\"],[LR[p],psi[p],psi[p].Conj(),M])\n",
    "                    \n",
    "                psi[p].set_labels(lbl1)\n",
    "                psi[p+1].set_labels(lbl2)\n",
    "                # Compute and save the environment tensor on the current site\n",
    "                LR[p+1] = anets[lr].Launch()\n",
    "                \n",
    "            print ('\\t',verbose[lr],'energy =',en)\n",
    "    return psi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90eb1049-bb27-4614-a35a-0398767e3e6f",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4.2 Find the ground state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d90dcb36-d452-4b4e-bbf6-c6768111348c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sweep 0 , chi=8\n",
      "\t [r->l] energy = -4.677035882069475\n",
      "\t [l->r] energy = -4.804469213722561\n",
      "Sweep 1 , chi=16\n",
      "\t [r->l] energy = -4.805131210889346\n",
      "\t [l->r] energy = -4.805151111016484\n",
      "Sweep 2 , chi=32\n",
      "\t [r->l] energy = -4.805151592051123\n",
      "\t [l->r] energy = -4.80515164604801\n",
      "Sweep 3 , chi=64\n",
      "\t [r->l] energy = -4.805151652062927\n",
      "\t [l->r] energy = -4.805151652838495\n",
      "Sweep 4 , chi=128\n",
      "\t [r->l] energy = -4.805151652926011\n",
      "\t [l->r] energy = -4.805151652935135\n",
      "Ground state energy for hz=0.1 is ~ -4.805151652932\n"
     ]
    }
   ],
   "source": [
    "Nsites = 20       # Number of sites\n",
    "\n",
    "# Get MPO tensor\n",
    "hz = 0.1          # Transverse field\n",
    "hx = 1e-8            # Longitudinal field\n",
    "M = TFIM_MPO_tensor (hz, hx)       # Make MPO tensor\n",
    "\n",
    "# Get initial state MPS\n",
    "psi = get_init_MPS (Nsites)\n",
    "\n",
    "# Define the bond dimensions for the sweeps\n",
    "maxdims = [8,16,32,64,128]\n",
    "cutoff = 1e-12\n",
    "\n",
    "# Run dmrg!!\n",
    "psi = dmrg (psi, M, maxdims, cutoff)\n",
    "\n",
    "print ('Ground state energy for hz=0.1 is ~ -4.805151652932')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fa00571-80cf-4965-9568-8760102c2c86",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 5. [Optional] Measure $\\langle S_z\\rangle$, $\\langle S_x\\rangle$, and the entanglement entropy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69a64a64-4f30-4e35-9236-cd5490466493",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 5.1 Define the measurement functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e083416-2735-4f79-9b7d-c6ad3c25f8fc",
   "metadata": {},
   "source": [
    "<img src=\"figures/szmea.png\" alt=\"drawing\" width=\"400\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd0ebf8a-5b9a-4c19-bc81-460bf884e18a",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Entanglement entropy\n",
    "#### $$ S_A = -\\mathrm{Tr}\\left(\\hat{\\rho}_A \\log \\hat{\\rho}_A\\right),$$\n",
    "#### where $\\hat{\\rho}_A$ is the reduced density matrix of subsystem A."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0812766c-aae5-419b-a673-d8e4e2a7ab3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def measure (MPS):\n",
    "    def mea_op (A, op):\n",
    "        anet = cytnx.Network()\n",
    "        anet.FromString([\"A: 1,2,3\",\n",
    "                         \"Op: 2,4\",\n",
    "                         \"A_Conj: 1,4,3\",\n",
    "                         \"TOUT:;\"])\n",
    "        anet.PutUniTensors ([\"Op\",\"A\",\"A_Conj\"],[op, A, A.Conj()])\n",
    "        return anet.Launch().item()\n",
    "\n",
    "    N = len(MPS)\n",
    "    Sz_ = cytnx.UniTensor (Sz, rowrank=1)\n",
    "    Sx_ = cytnx.UniTensor (Sx, rowrank=1)                                                                                                                                                                                                                                                                                                                                                                                                                 \n",
    "\n",
    "    psi = MPS.copy()\n",
    "    Ss = np.full (N-1, np.nan)\n",
    "    szs = np.full (N, np.nan)\n",
    "    sxs = np.full (N, np.nan)\n",
    "\n",
    "    for p in range(Nsites-2,-1,-1):\n",
    "        psi[p+1].set_rowrank(1)\n",
    "        # Measure Sz\n",
    "        szs[p+1] = mea_op (psi[p+1], Sz_)\n",
    "        sxs[p+1] = mea_op (psi[p+1], Sx_)\n",
    "\n",
    "        # Move the orthogonality center to the left site\n",
    "        s, u, vt = cytnx.linalg.Svd (psi[p+1])\n",
    "        psi[p].put_block(cytnx.Contract(psi[p], cytnx.Contract(u, s)).get_block())\n",
    "        psi[p+1].put_block(vt.get_block())\n",
    "        \n",
    "        # Measure the entanglement entropy\n",
    "        S = 0\n",
    "        for si in s.get_block_():\n",
    "            rhoi = si.item()**2\n",
    "            S += -rhoi * np.log(rhoi)\n",
    "        Ss[p] = S\n",
    "    psi[0].set_rowrank(1)\n",
    "    szs[0] = mea_op (psi[0], Sz_)\n",
    "    sxs[0] = mea_op (psi[0], Sx_)\n",
    "\n",
    "    return szs, sxs, Ss"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24d4b4d6-58d7-4f43-81b0-b4e7ddd5efe4",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 5.2 Do the measurements and save files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "332d88e5-ef0f-4d90-852c-1ba251ec39f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "szs, sxs, Ss = measure (psi)\n",
    "np.save(\"datas/szs\",szs)\n",
    "np.save(\"datas/sxs\",sxs)\n",
    "np.save(\"datas/Ss\",Ss)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec706bcb-e5d4-4022-922f-2891210682a7",
   "metadata": {},
   "source": [
    "## Homework 1.\n",
    "In the example above we've study the traverse field Ising model at small $h_z$, and find $\\langle S_z \\rangle = 0 $ and $\\langle S_x \\rangle = 0.5 $, how about at **large** $h_z$ like $h_z = 2$ ? Can you find the **transition point** ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f56de58-4c14-4d0b-bc7e-86ffa317e728",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "294f04a4-362f-453d-a0fc-1b36458c9f55",
   "metadata": {},
   "source": [
    "## Homework 2.\n",
    "Find the ground state energy of the **XX model** given by the Hamiltonian:\n",
    "\n",
    "$$\\hat{H} = \\sum_{i=1}^{N-1} -\\hat{S}^x_i\\hat{S}^x_{i+1} -  \\hat{S}^y_i\\hat{S}^y_{i+1}$$\n",
    "\n",
    "You can use the following function to get the exact value given the desired system size $N$, use it to **benchmark** your DRMG result with $N=16, 32,100$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "04fcad05-0ecc-4a84-9338-038abe3496d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def getXXexact(N):\n",
    "    import numpy as np\n",
    "    H = np.diag(np.ones(N-1),k=1) + np.diag(np.ones(N-1),k=-1)\n",
    "    D = LA.eigvalsh(H)\n",
    "    EnExact = 2*sum(D[D < 0])\n",
    "    return EnExact\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13dd9bd7-d9b0-497b-a3c4-67e970cffa9b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "91dba0e0-f2b3-4181-a3a6-46317834f999",
   "metadata": {},
   "source": [
    "## Homework 3.\n",
    "Find the ground state energy of the **antiferromagnet Heisenberg model** given by the Hamiltonian:\n",
    "\n",
    "$$\\hat{H} = \\sum_{i=1}^{N-1} \\hat{S}^x_i\\hat{S}^x_{i+1} + \\hat{S}^y_i\\hat{S}^y_{i+1} + \\hat{S}^z_i\\hat{S}^z_{i+1}$$\n",
    "\n",
    "Please try several values of $N$ and **extraplolated** to $N \\to \\infty$ and compare to the exact value at **thermaldynamic limit**:\n",
    "$$1/4 - log(2) \\approx -0.4431$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "267531c0-fd2e-4aa8-b469-084c764b5910",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}