2 - Simulation Packages Tour.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we'll go through some examples of quantum physics simulation packages in Julia. We'll look at tensor network simulation packages, eigensolving packages, and leave with some ideas for projects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tensor network simulation\n",
    "\n",
    "Tensor network methods are a family of techniques commonly employed in many body physics simulation. We saw an example already at this school in DMRG, but there are of course many others (and many types of DMRG), such as PEPS methods, tree tensor networks, MERA, etc.\n",
    "\n",
    "ITensor is a package that was originally written in C++ but was recently ported from the ground up to Julia. We'll go through some of its features that make it a good tool for physics simulation. Note that there are other tensor network simulation packages in Julia! Jutho Hagaeman has a nice selection on [his GitHub](https://github.com/jutho/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basics: building some tensors, doing operations to them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first, let's install the package!\n",
    "using Pkg; Pkg.add(\"ITensors\")\n",
    "using ITensors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ITensors have indices, which have descriptive \"tags\"\n",
    "a = Index(2, \"a\")\n",
    "b = Index(3, \"b\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can build up an ITensor from the indices\n",
    "A = ITensor(a, b)\n",
    "# and set its elements\n",
    "A[a=>1, b=>3]=2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the default element type of an ITensor is Float64\n",
    "@show eltype(A)\n",
    "# but we can control this\n",
    "B = ITensor(ComplexF64, b, a)\n",
    "@show eltype(B)\n",
    "B[a=>1, b=>3] = 2*im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can contract two tensors using the * operator\n",
    "A * B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ITensors.jl will automatically figure out the output indices for us:\n",
    "c = Index(2, \"c\")\n",
    "d = Index(3, \"d\")\n",
    "C = randomITensor(a, b, c)\n",
    "D = randomITensor(c, d, a)\n",
    "@show inds(C*D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ITensors also provides decomposition methods like svd and eigen\n",
    "svd(C, a) # a here is treated as a \"left\"/\"row\" index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svd(C, a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building MPS and MPO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ITensor provides methods to automatically generate states with\n",
    "# the correct index flux and dimensionality for various problems\n",
    "N     = 10\n",
    "sites = siteinds(\"S=1/2\",N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# and we can build MPS from these\n",
    "ψ     = randomMPS(sites,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# or MPO\n",
    "O     = randomMPO(sites,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ITensors will automatically handle matching up the correct \"Site\" indices\n",
    "O * ψ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can easily compute inner products\n",
    "ϕ     = randomMPS(sites,10)\n",
    "dot(ϕ, O, ψ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are quite a few other MPS/MPO methods you can check out:\n",
    "- `inner`\n",
    "- `projMPS`\n",
    "- `projMPO`\n",
    "- Try constructing with different `siteinds` (e.g. `S=1`, fermions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DMRG\n",
    "We don't want to have to construct the MPO element by element -- how annoying! Instead, `ITensors` provides a way to automatically construct MPOs from equations -- `AutoMPO`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10\n",
    "ampo = AutoMPO()\n",
    "h = 0.2\n",
    "for j=1:N-1\n",
    "    add!(ampo, -1,\"Sz\",j,\"Sz\",j+1)\n",
    "end\n",
    "for j=1:N\n",
    "    add!(ampo, h,\"Sx\",j)\n",
    "end\n",
    "# Convert these terms to an MPO tensor network\n",
    "sites = siteinds(\"S=1/2\",N)\n",
    "H = MPO(ampo,sites)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# With this MPO, we can do DMRG\n",
    "ψ₀ = randomMPS(sites,10)\n",
    "sweeps = Sweeps(15)\n",
    "setmaxdim!(sweeps, 10,20,100,100,200) # set the maximum bond dimension at each sweep\n",
    "setcutoff!(sweeps, 1E-10) # set truncation error cutoff (for all sweeps)\n",
    "energy, ψ = dmrg(H,ψ₀,sweeps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also build \"observers\" into our DMRG to track convergence over time or observables of interest\n",
    "Sz_observer = DMRGObserver([\"Sz\"],sites,energy_tol=1E-7)\n",
    "energy, ψ = dmrg(H,ψ₀,sweeps, observer = Sz_observer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take some time to explore. Some suggestions:\n",
    "- Try the [examples](https://github.com/ITensor/ITensors.jl/tree/master/examples/dmrg) in the ITensors repo (you might find the Hubbard models particularly interesting, or models with fermions).\n",
    "- Play around with some of the [options](https://itensor.github.io/ITensors.jl/dev/DMRG.html) and varities of DMRG that ITensors supports\n",
    "- Try changing up the `AutoMPO` to include disorder. You might find [Distributions.jl](https://github.com/JuliaStats/Distributions.jl) helpful for generating certain distributions of random numbers.\n",
    "- Try the [multithreading](https://itensor.github.io/ITensors.jl/dev/Multithreading.html) support. Do you notice a speed difference?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Krylov methods\n",
    "\n",
    "Krylov methods are used to perform eigensolving or time evolution. You likely met them already in the lecture on exact diagonalization. We'll use [KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl) to perform some basic eigensolving. Note again that there are [*quite a few*](https://jutho.github.io/KrylovKit.jl/latest/#Package-features-and-alternatives) Julia packages aiming to solve these types of problems, and you're encouraged to check them out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first, let's install the package!\n",
    "using Pkg; Pkg.add(\"KrylovKit\")\n",
    "using KrylovKit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# KrylovKit provides methods for (extremal) eigenvalue solving, operator exponentiation, SVD, and several others.\n",
    "# Let's try a simple example with an explicit matrix\n",
    "n = 1024 # small\n",
    "A = rand(ComplexF64,(n,n)) .- one(ComplexF64)/2\n",
    "A = (A+A')/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ishermitian(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tol = 10*n*eps(real(ComplexF64))\n",
    "λ, ϕ, info = eigsolve(A, 1, :SR, orth=ClassicalGramSchmidt(), tol=tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# KrylovKit also supports \"matrix-free methods\"\n",
    "# for these we need to provide a starting guess x₀\n",
    "# Here, I'm going to be lazy and provide a \"fake\" method\n",
    "f(v::AbstractVector) = A*v\n",
    "x₀ = randn(n)\n",
    "tol = 10*n*eps(real(ComplexF64))\n",
    "λ, ϕ, info = eigsolve(f, x₀, 1, :SR, orth=ClassicalGramSchmidt(), tol=tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's try the exponentiator\n",
    "x = randn(n)\n",
    "alg = Arnoldi(orth = ClassicalGramSchmidt(), krylovdim = n, maxiter = 2, tol = 10*n*eps(real(ComplexF64)))\n",
    "y = exponentiate(A, 1, x, alg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = exponentiate(A, 2im/π, x, alg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = exponentiate(f, 2im/π, x, alg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take some time to explore. Some suggestions:\n",
    "- Try different solver methods -- SVD, for example\n",
    "- Try building your own matrices or matrix-free functions using what you learned in David Luitz's lecture\n",
    "- How do different orthogonalizers affect things?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.1-pre",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}