diff --git a/exercises/ex13_skeleton/heisenberg_dmrg_skeleton.ipynb b/exercises/ex13_skeleton/heisenberg_dmrg_skeleton.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..88231f767a51df5913c7e8c26d55bb450d403343 --- /dev/null +++ b/exercises/ex13_skeleton/heisenberg_dmrg_skeleton.ipynb @@ -0,0 +1,222 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "from numpy import transpose as tr, conjugate as co\n", + "from scipy.linalg import expm, svd\n", + "from scipy.sparse.linalg import eigsh, LinearOperator\n", + "import math\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "plt.rcParams['figure.figsize'] = 16, 9" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Useful functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def dot(A,B):\n", + " \"\"\" Does the dot product like np.dot, but preserves the shapes also for singleton dimenstions \"\"\"\n", + " s1 = A.shape\n", + " s2 = B.shape\n", + " return np.dot(A,B).reshape((s1[0],s2[1]))\n", + "\n", + "def left_canonize(M1,M2,return_S = False):\n", + " \"\"\" Left normalizes M1 into A matrix, M2 looses its canonization\"\"\"\n", + " s, da, db = M1.shape\n", + " U, S, Vh = svd(M1.reshape((s*da,db)))\n", + " U = U.reshape((s,da,U.shape[1]))[:,:,:db] # is now A matrix\n", + " M2 = np.tensordot(dot(np.diag(S[:db]),Vh[:db,:]),M2,axes=((1),(1))) # looses its canonization\n", + " if return_S:\n", + " return U, M2, S\n", + " else:\n", + " return U, M2\n", + "\n", + "def right_canonize(M1,M2,return_S = False):\n", + " \"\"\" Right normalizes M2 into B matrix, M1 looses its canonization \"\"\"\n", + " s, da, db = M2.shape\n", + " U, S, Vh = svd(M2.transpose((1,0,2)).reshape((da,s*db)))\n", + " Vh = Vh.reshape((Vh.shape[0],s,db)).transpose((1,0,2))[:,:da,:] # is now B matrix\n", + " M1 = np.tensordot(M1,dot(U[:,:da],np.diag(S[:da])),axes=((2),(0))) # looses its canonization\n", + " if return_S:\n", + " return M1, Vh, S\n", + " else:\n", + " return M1, Vh\n", + " \n", + "def add_site_to_R(Rexp,M,W):\n", + " ## to be completed\n", + " return Rexp\n", + "\n", + "def add_site_to_L(Lexp,M,W):\n", + " ## to be completed\n", + " return Lexp\n", + "\n", + "def Hpsi(L,W,R,M):\n", + " \"\"\"\n", + " Calculates the \"matrix-vector product\" at arbitrary site from \n", + " L ... L-expression, belonging to the \"matrix\"\n", + " W ... W-array at site l, belonging to the \"matrix\"\n", + " R ... R-expression, belonging to the \"matrix\"\n", + " M ... M-array, the \"vector\"\n", + " \"\"\"\n", + " ## to be completed" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Definitions and construction of MPO" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "J = 1.\n", + "L = 40 # Length of chain\n", + "s = 2 # Local dimension of Hilbert space\n", + "chi = 60 # The maximum matrix dimension, from which on the matrices are truncated\n", + "nmax = 10000 # Maximum number of iterations\n", + "\n", + "## Construct the MPO for the Heisenberg model below\n", + "# ..." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Random guess as groundstate" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "## random starting configuration\n", + "M = [np.random.standard_normal((s,min(chi,s**i,s**(L-i)),min(chi,s**(i+1),s**(L-1-i)))) for i in range(L)]\n", + "\n", + "## need a right canonized MPS\n", + "for i in range(L-1,0,-1):\n", + " M[i-1], M[i] = right_canonize(M[i-1],M[i])\n", + "## normalization, throw away the left part\n", + "_, M[0] = right_canonize(np.ones((1,1,1)),M[0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### DMRG sweeps" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "## get initial R expression\n", + "Rexp = [np.ones((1,1,1))]\n", + "for i in range(1,L):\n", + " Rexp.append(add_site_to_R(Rexp[i-1],M[-i],W[-i]))\n", + " \n", + "enold = 0\n", + " \n", + "for n in range(nmax):\n", + " ## right sweep\n", + " Lexp = [np.ones((1,1,1))]\n", + " for i in range(L-1):\n", + " print('.', end=\"\")\n", + " # ----define a linear operator with square dimension s*D^2 via providing a matrix vector product---\n", + " dD2 = s*Lexp[i].shape[0]*Rexp[L-1-i].shape[0]\n", + " Hop = LinearOperator((dD2,dD2),matvec= lambda M: Hpsi(Lexp[i],W[i],Rexp[L-1-i],M))\n", + " # ----Use an iterative eigenvalue solver, M[i] needs to be flattened into a vector---\n", + " en,V = eigsh(Hop,k=1,v0=M[i].flatten(),tol=1e-3 if n < 2 else 0)\n", + " # ----Reshape M[i] back into its 3D array shape----\n", + " M[i] = V.reshape((s,Lexp[i].shape[0],Rexp[L-1-i].shape[0]))\n", + " # ----Left normalize M[i]---\n", + " M[i], M[i+1] = left_canonize(M[i],M[i+1])\n", + " # ----add site to L-expression----\n", + " Lexp.append(add_site_to_L(Lexp[i],M[i],W[i]))\n", + " print('')\n", + "\n", + " ## left sweep\n", + " Rexp = [np.ones((1,1,1))]\n", + " for i in range(L-1,0,-1):\n", + " print('.', end=\"\")\n", + " ## to be completed\n", + " print('')\n", + " \n", + " # === check convergence ===\n", + " print(n,'dE = ',abs(en-enold))\n", + " if abs(en-enold) < 1e-8:\n", + " print(\"Converged after \"+str(n)+\" sweeps!\")\n", + " print(\"Ground-state energy: \" + str(en))\n", + " break\n", + " enold = en\n", + " # === check convergence ===\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "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.5.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/exercises/ex13_skeleton/heisenberg_imaginary_time_evolution.ipynb b/exercises/ex13_skeleton/heisenberg_imaginary_time_evolution.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..47a29f72da6cecd3d7f8dce5f2cf22fa9b3645a2 --- /dev/null +++ b/exercises/ex13_skeleton/heisenberg_imaginary_time_evolution.ipynb @@ -0,0 +1,364 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "from numpy import transpose as tr, conjugate as co\n", + "from scipy.linalg import expm, svd\n", + "import math\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "plt.rcParams['figure.figsize'] = 16, 9" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Some useful functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def dot(A,B):\n", + " \"\"\" Does the dot product like np.dot, but preserves the shapes also for singleton dimenstions \"\"\"\n", + " s1 = A.shape\n", + " s2 = B.shape\n", + " return np.dot(A,B).reshape((s1[0],s2[1]))\n", + "\n", + "def truncated_svd(thetas,chi):\n", + " \"\"\" Does an svd on two-site matrix thetas, and truncates to chi or the last nonzero singular value \"\"\"\n", + " U, S, Vh = svd(thetas,full_matrices = False)\n", + " # trunkieren\n", + " ind = np.where(np.isclose(np.cumsum(S[::-1])[::-1],0))[0]\n", + " if len(ind)>0:\n", + " chi_tr = min(chi,max(1,ind[0]))\n", + " else:\n", + " chi_tr = chi\n", + " S=S[:chi_tr]\n", + " w=1.-np.sum(S**2)\n", + " S/=math.sqrt(1-w)\n", + " U=U[:,:chi_tr]\n", + " Vh=Vh[:chi_tr,:]\n", + " return U, S, Vh\n", + "\n", + "def left_canonize(thetas,s,chi, return_S = False):\n", + " \"\"\" Splits up a two-site matrix thetas into two one-site matrices such that the left one is canonized \"\"\"\n", + " da, dg = thetas.shape[2], thetas.shape[3]\n", + " thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes\n", + " U, S, Vh = truncated_svd(thetas,chi)\n", + "\n", + " db = len(S)\n", + " U = U.reshape((da,s,db)).transpose((1,0,2))\n", + " Vh = Vh.reshape((db,dg,s)).transpose((2,0,1)) \n", + " for s2 in range(s):\n", + " Vh[s2] = dot(np.diag(S),Vh[s2])\n", + " if return_S:\n", + " return U, Vh, S\n", + " else:\n", + " return U, Vh\n", + "\n", + "def right_canonize(thetas,s,chi, return_S = False):\n", + " \"\"\" Splits up a two-site matrix thetas into two one-site matrices such that the right one is canonized \"\"\"\n", + " da, dg = thetas.shape[2], thetas.shape[3]\n", + " thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes\n", + " U, S, Vh = truncated_svd(thetas,chi)\n", + "\n", + " db = len(S)\n", + " U = U.reshape((da,s,db)).transpose((1,0,2))\n", + " Vh = Vh.reshape((db,dg,s)).transpose((2,0,1)) \n", + " for s1 in range(s):\n", + " U[s1] = dot(U[s1],np.diag(S))\n", + " if return_S:\n", + " return U, Vh, S\n", + " else:\n", + " return U, Vh\n", + "\n", + "\n", + "def apply_two_site_H(theta,H,s):\n", + " \"\"\" Applies a two-site operator on the two-site matrix theta \"\"\"\n", + " thetas = np.zeros_like(theta)\n", + " for s1p in range(s):\n", + " for s2p in range(s):\n", + " for s1 in range(s):\n", + " for s2 in range(s):\n", + " thetas[s1p,s2p] += H[s1p*s+s2p,s1*s+s2] * theta[s1,s2]\n", + "\n", + " return thetas\n", + "\n", + "def combine_two_matrices(M1,M2,s):\n", + " \"\"\" Combines the two neighbouring one-site matrices M1 and M2 into a two-site matrix theta \"\"\"\n", + " da, db, dg=M1.shape[1], M1.shape[2], M2.shape[2]\n", + " assert M1.shape[2] == M2.shape[1]\n", + "\n", + " theta = np.zeros((s,s,da,dg),dtype=complex)\n", + " for s1 in range(s):\n", + " for s2 in range(s):\n", + " theta[s1,s2]=dot(M1[s1],M2[s2])\n", + "\n", + " return theta" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Definition of system and initializiation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "J = 1.\n", + "L = 20 # Length of chain\n", + "s = 2 # Local dimension of Hilbert space\n", + "dt = 0.05 # The imaginary timestep (should be 0.01 or smaller for good accuracy)\n", + "chi = 60 # The maximum matrix dimension, from which on the matrices are truncated\n", + "nmax = 10000 # Maximum number of iterations\n", + "nskip = 10 # Check energy convergence after nskip steps\n", + "\n", + "# Two-site Hamiltonian\n", + "H = np.array([[J/4,0,0,0],\n", + " [0,-J/4,J/2,0],\n", + " [0,J/2,-J/4,0],\n", + " [0,0,0,J/4]])\n", + "\n", + "# Imaginary time evolution operator\n", + "exp_beta_H=expm(-H*dt);\n", + "\n", + "# antiferromagnetic starting configuration\n", + "M = []\n", + "for i in range(L):\n", + " ar = np.zeros((2,1,1),dtype=complex)\n", + " ar[0,0,0] = i%2\n", + " ar[1,0,0] = (i+1)%2\n", + " M.append(ar)\n", + "\n", + "even = np.array(range(0,L-1,2))\n", + "odd = np.array(range(L+L%2-2,0,-2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Imaginary time evolution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "enold = 0. # for checking energy convergence\n", + "\n", + "for n in range(nmax):\n", + " # ++++ Trotter scheme: first even bonds, then odd bonds\n", + " # === apply exp_beta_H on all even bonds ===\n", + " for j in even:\n", + " # go from left to right\n", + " theta = combine_two_matrices(M[j],M[j+1],s)\n", + " thetas = apply_two_site_H(theta,exp_beta_H,s) \n", + " M[j], M[j+1] = left_canonize(thetas,s,chi)\n", + "\n", + " # advance left canonization by a further step\n", + " if j < L-2:\n", + " theta = combine_two_matrices(M[j+1],M[j+2],s) \n", + " M[j+1], M[j+2] = left_canonize(theta,s,chi)\n", + " # === apply exp_beta_H on all even bonds ===\n", + "\n", + " # === renormalize the state on the last site ===\n", + " da = M[L-1].shape[1]\n", + " theta = M[L-1].reshape((s*da,1))\n", + " U, _, _ = truncated_svd(theta,chi) # throw away norm\n", + " M[L-1] = U.reshape((s,da,1))\n", + " # === renormalize the state on the last site ===\n", + "\n", + " if L%2 == 0:\n", + " # the right-most matrix has to be right canonized in this case\n", + " theta = combine_two_matrices(M[L-2],M[L-1],s)\n", + " M[L-2], M[L-1] = right_canonize(theta,s,chi)\n", + "\n", + " # === apply exp_beta_H on all odd bonds ===\n", + " for j in odd:\n", + " # go from right to left\n", + " theta = combine_two_matrices(M[j-1],M[j],s)\n", + " thetas = apply_two_site_H(theta,exp_beta_H,s)\n", + " M[j-1], M[j] = right_canonize(thetas,s,chi)\n", + "\n", + " # advance right canonization by a further step\n", + " if j>1:\n", + " theta = combine_two_matrices(M[j-2],M[j-1],s)\n", + " M[j-2], M[j-1] = right_canonize(theta,s,chi)\n", + " # === apply exp_beta_H on all odd bonds ===\n", + "\n", + " # === renormalize the state on the first site ===\n", + " db = M[0].shape[2]\n", + " theta = M[0].reshape((1,s*db))\n", + " _, _, Vh = truncated_svd(theta,chi) # throw away norm\n", + " M[0] = Vh.reshape((s,1,db))\n", + " # === renormalize the state on the first site ===\n", + "\n", + " # ++++ Calculate energy and check convergence ++++\n", + " if n%nskip == 0:\n", + " # Measure Energy via sum of two-site Operators H\n", + " en = 0.\n", + " for j in range(L-1): \n", + " theta = combine_two_matrices(M[j],M[j+1],s)\n", + " thetas = apply_two_site_H(theta,H,s) \n", + " for s1 in range(s):\n", + " for s2 in range(s):\n", + " en += np.trace(dot(thetas[s1,s2],co(tr(theta[s1,s2])))).real\n", + " M[j], M[j+1] = left_canonize(theta,s,chi) # Proceed with left canonization\n", + "\n", + " # make right canonized MPS for next tDMRG step\n", + " for j in range(L-1,0,-1):\n", + " theta = combine_two_matrices(M[j-1],M[j],s)\n", + " M[j-1], M[j] = right_canonize(theta,s,chi)\n", + "\n", + " print(n,'dE = ',abs(en-enold))\n", + " if abs(en-enold) < 1e-8:\n", + " print(\"Converged after \"+str(n)+\" steps!\")\n", + " print(\"Ground-state energy: \" + str(en))\n", + " break\n", + " enold = en\n", + " # === check convergence ===\n", + "\n", + "if n == nmax:\n", + " print(\"Maximum iterations reached! Convergence criterium not satisfied!\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Implement the Heisenberg MPO below\n", + "\n", + "You can build it as a list of four-dimensional arrays with indexing $\\hat W^{[i]} = W_{\\sigma'_i\\, \\sigma_i\\, b_{l-1}\\, b_l}$. It can be represented as \"matrices\" with indices $(b_{l-1}, b_l)$ which have single-site operators as entries having $(\\sigma'_l,\\sigma_l)$ indices. For the Heiseberg model \n", + " \\begin{equation}\n", + " H = \\sum_{i=1}^{L-1} \\left( \\frac{J}{2} S_i^+ S_{i+1}^- + \\frac{J}{2} S_i^- S_{i+1}^+ + J S_i^z S_{i+1}^z\\right) ,\n", + " \\end{equation}\n", + " the matrices are\n", + " \\begin{equation}\n", + " \\hat W^{[i]} = \\begin{bmatrix}\n", + " \\hat I & 0 & 0 & 0 & 0 \\\\\n", + " S^+_i & 0 & 0 & 0 & 0 \\\\\n", + " S^-_i & 0 & 0 & 0 & 0 \\\\\n", + " S^z_i & 0 & 0 & 0 & 0 \\\\\n", + " 0 & (J/2) S^-_i & (J/2) S^+_i & J S^z_i & \\hat I\n", + " \\end{bmatrix},\n", + " \\end{equation}\n", + " and for the first and last site\n", + " \\begin{equation}\n", + " \\hat W^{[1]} = \\begin{bmatrix}\n", + " 0 & (J/2)S^-_0 & (J/2) S^+_0 & J S^z_0 & \\hat I\n", + " \\end{bmatrix}, \\ \\ \\hat W^{[L]} = \\begin{bmatrix}\n", + " \\hat I \\\\ S^+_L \\\\ S^-_L \\\\ S^z_L \\\\ 0\n", + " \\end{bmatrix}.\n", + " \\end{equation}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "## Generate MPO\n", + "W = []\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "### Build iteratively the L- or R-expression and evaluate the energy of the MPS using above MPO\n", + "\n", + "You can use a three-dimensional array for $R$ with the indexing $R_{a'_l\\, b_l\\, a_l}$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def add_site_to_R(Rexp,M,W):\n", + " \"\"\"Rexp... 3-d array\n", + " M ... M-array of the site to be added\n", + " W ... W-array of the site to be added\"\"\"\n", + " ## to be completed\n", + " return Rexp\n", + "\n", + "def add_site_to_L(Lexp,M,W):\n", + " ## to be completed\n", + " return Lexp\n", + "\n", + "## Calculate energy via iteratively constructing the R-expression\n", + "Rexp = np.ones((1,1,1))\n", + "for i in range(1,L+1):\n", + " Rexp = add_site_to_R(Rexp,M[-i],W[-i])\n", + "print('Energy from R-expression: ',np.squeeze(Rexp)) # should be total energy\n", + "\n", + "## Calculate energy via iteratively constructing the L-expression\n", + "# ..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "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.5.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/exercises/exercise13.pdf b/exercises/exercise13.pdf new file mode 100644 index 0000000000000000000000000000000000000000..26c4a278490f80cf548d37f35dc8fd3880248f12 Binary files /dev/null and b/exercises/exercise13.pdf differ