DMRG_Ising.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Linear Algebra\n",
    "import numpy as np\n",
    "from scipy.linalg import svd\n",
    "from numpy.linalg import qr\n",
    "import scipy.sparse.linalg.eigen.arpack as arp\n",
    "\n",
    "from timeit import timeit\n",
    "\n",
    "# Python Plotting'\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Own modules\n",
    "from DMRG_022_LP import * # here is the main part of the DMRG code\n",
    "\n",
    "# we'll need them a lot\n",
    "sxx = np.array([[0., 1.], [1., 0.]])\n",
    "syy = np.array([[0., -1j], [1j, 0.]])\n",
    "szz = np.array([[1., 0.], [0., -1.]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class M_MPO(object):\n",
    "    '''\n",
    "    MPO for the total magnetization operator M = \\\\sum_i \\\\sigma_i^z\n",
    "    '''\n",
    "    def __init__(self,L):\n",
    "        self.L = L\n",
    "        self.d = 2\n",
    "\n",
    "        self.sz = np.array([[1., 0.], [0., -1.]])\n",
    "        self.sx = np.array([[0.,1.],[1.,0.]])\n",
    "        self.id = np.eye(2)\n",
    "        self.Ws = self.init_W()  # [?] in toycode they do this a bit strange\n",
    "                                #and make the defintion self.W inside init_W, is this necessary or rather 'style?'\n",
    "    \n",
    "    def init_W(self):\n",
    "        # Create the matrix of operators (2x2 matrices)\n",
    "        # (virutal left, virtual right, physical out, physical in)\n",
    "        # vL, vR, s, s'\n",
    "        Ws = []\n",
    "        # First vector\n",
    "        w = np.zeros((1,2,self.d,self.d),dtype=np.complex)\n",
    "        w[0,0] = self.id\n",
    "        w[0,1] = self.sx\n",
    "        Ws.append(w)\n",
    "        # Ws[1:L-2]\n",
    "        w = np.zeros((2,2,self.d,self.d),dtype=np.complex)\n",
    "        w[0,0] = self.id\n",
    "        w[0,1] = self.sx\n",
    "        w[1,1] = self.id\n",
    "        # create L-2 times the same matrix\n",
    "        for i in range(1,self.L-1):\n",
    "            Ws.append(w)\n",
    "        # Last vector\n",
    "        w = np.zeros((2,1,self.d,self.d),dtype=np.complex)\n",
    "        w[0,0] = self.sx\n",
    "        w[1,0] = self.id\n",
    "        Ws.append(w)\n",
    "        return Ws\n",
    "\n",
    "def dmrg_m_MPO(L,J,g,chi_max=50,conf=1e-4,test=False,eps=-1,tol=0,maxiter=None):\n",
    "    '''\n",
    "    calculate mean magnetization m = <\\\\sum_i \\\\sigma_i^z>\n",
    "    FROM MPO\n",
    "    '''\n",
    "    # Perform DMRG\n",
    "    ## Initialize random MPS and MPO\n",
    "    init_mps = random_MPS(L=L,d=2,chi_max=chi_max)\n",
    "    mpo = Ising_MPO(L=L,d=2,h=g,J=J)\n",
    "    dmrg = DMRG(init_mps,mpo,eps=eps,chi_max=chi_max,test=False,tol=tol,maxiter=maxiter)\n",
    "    ## Left- and Right-sweep\n",
    "    diff = 5.*conf # just to make it bigger than conf\n",
    "    counter = 0\n",
    "    while diff > conf:\n",
    "        if test:\n",
    "            print(counter,diff)\n",
    "        counter += 1\n",
    "        ### note: we could also simply use the current dmrg.e value to speed things up. This is just to give justification to the expect_mpo(mps,mpo) function\n",
    "        dmrg.left_to_right() # appends the current e-values for each site\n",
    "        e_in = dmrg.e[-1]\n",
    "        dmrg.right_to_left()\n",
    "        e_end = dmrg.e[-1]\n",
    "        diff = abs(e_in - e_end)/abs(e_in)\n",
    "        #if test and counter > 1:\n",
    "        #    print('diff = ',diff)\n",
    "        if counter > 5:\n",
    "            print('Exceeded 5 epochs')\n",
    "            break\n",
    "    if test:\n",
    "        print('g = ',g,' - epochs = ',counter-1,' - rel. diff = ',100*diff,'%')\n",
    "        print('E_dmrg = ',dmrg.e[-1])\n",
    "    # Calculate magnetization\n",
    "    return expect_mpo(dmrg.MPS,M_MPO(L))/L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 4.9999999999999996e-06\n",
      "g =  -0.0  - epochs =  0  - rel. diff =  8.943037881118471e-13 %\n",
      "E_dmrg =  -28.999999999999844\n",
      "0 4.9999999999999996e-06\n",
      "g =  -0.2  - epochs =  0  - rel. diff =  6.058338298401714e-14 %\n",
      "E_dmrg =  -29.320859151574325\n",
      "0 4.9999999999999996e-06\n",
      "g =  -0.4  - epochs =  0  - rel. diff =  1.1727365734691711e-14 %\n",
      "E_dmrg =  -30.294217466851215\n",
      "0 4.9999999999999996e-06\n",
      "g =  -0.6000000000000001  - epochs =  0  - rel. diff =  4.427465709733018e-07 %\n",
      "E_dmrg =  -31.95661630391031\n",
      "0 4.9999999999999996e-06\n"
     ]
    }
   ],
   "source": [
    "gs = np.linspace(0.,2.,11)\n",
    "mag_MPO_30 = np.array([dmrg_m_MPO(L=30,J=-1.,g=-1*g,chi_max=40,test=True,conf=1e-6,eps=-1) for g in gs])\n",
    "mag_MPO_50 = np.array([dmrg_m_MPO(L=50,J=-1.,g=-1*g,chi_max=40,test=True,conf=1e-6,eps=-1) for g in gs])\n",
    "mag_MPO_70 = np.array([dmrg_m_MPO(L=70,J=-1.,g=-1*g,chi_max=40,test=True,conf=1e-6,eps=-1) for g in gs])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savez('data/dmrg_MPO_L.npz',mag_MPO_30 = mag_MPO_30,mag_MPO_50 = mag_MPO_50,mag_MPO_70 = mag_MPO_70)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data=np.load('data/dmrg_MPO_L.npz')\n",
    "mag_MPO_30,mag_MPO_50,mag_MPO_70 = data['mag_MPO_30'],data['mag_MPO_50'],data['mag_MPO_70']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=1,figsize=(10,6))\n",
    "### changed minus in MPO to check something\n",
    "ax=axes\n",
    "ax.plot(gs,mag_MPO_30,'-',label='L = 30')\n",
    "ax.plot(gs,mag_MPO_50,'-',label='L = 50')\n",
    "ax.plot(gs,mag_MPO_70,'-',label='L = 70')\n",
    "ax.set(xlabel='g', ylabel='m',\n",
    "       title='chi_max=40,conf=1e-6,eps=1e-3')\n",
    "ax.legend()\n",
    "#fig.savefig(\"mag_comp.png\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}