ed_sz_conserved_1d_heisenberg_jl.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "snoob (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function snoob(x)\n",
    "    next = 0\n",
    "    if(x>0)\n",
    "        smallest = x & -(x)\n",
    "        ripple = x + smallest\n",
    "        ones = xor(x,ripple)\n",
    "        ones = div((ones >> 2),smallest)\n",
    "        next = ripple | ones\n",
    "    end\n",
    "    return next\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count_bit (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function count_bit(n)\n",
    "    count = 0\n",
    "    while (n>0)\n",
    "        count += n & 1\n",
    "        n >>= 1\n",
    "    end\n",
    "    return count \n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "init_parameters (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function init_parameters(N,Sz)\n",
    "    Nup = div(N,2) + Sz\n",
    "    Nhilbert = binomial(N,Nup)\n",
    "    ihfbit = 1 << (div(N,2))\n",
    "    irght = ihfbit-1\n",
    "    ilft = xor(((1<<N)-1),irght)\n",
    "    iup = (1<<(N-Nup))-1\n",
    "    return Nup, Nhilbert, ihfbit, irght, ilft, iup\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "make_list (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function make_list(N,Nup,Nhilbert,ihfbit,irght,ilft,iup)\n",
    "    list_1 = zeros(Int,Nhilbert)\n",
    "    list_ja = zeros(Int,ihfbit)\n",
    "    list_jb = zeros(Int,ihfbit)\n",
    "    ii = iup\n",
    "    ja = 0\n",
    "    jb = 0\n",
    "    ia_old = ii & irght\n",
    "    ib_old = div((ii & ilft),ihfbit)\n",
    "    list_1[1] = ii\n",
    "    list_ja[ia_old+1] = ja\n",
    "    list_jb[ib_old+1] = jb\n",
    "    ii = snoob(ii)\n",
    "    for i in 1:Nhilbert-1\n",
    "        ia = ii & irght\n",
    "        ib = div((ii & ilft),ihfbit)\n",
    "        if ib == ib_old\n",
    "            ja += 1\n",
    "        else\n",
    "            jb += ja+1\n",
    "            ja = 0\n",
    "        end\n",
    "        list_1[i+1] = ii\n",
    "        list_ja[ia+1] = ja\n",
    "        list_jb[ib+1] = jb\n",
    "        ia_old = ia\n",
    "        ib_old = ib\n",
    "        ii = snoob(ii)\n",
    "    end\n",
    "    return list_1, list_ja, list_jb\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_ja_plus_jb (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function get_ja_plus_jb(ii,irght,ilft,ihfbit,list_ja,list_jb)\n",
    "    ia = ii & irght\n",
    "    ib = div((ii & ilft),ihfbit)\n",
    "    ja = list_ja[ia+1]\n",
    "    jb = list_jb[ib+1]\n",
    "    return ja+jb\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "make_hamiltonian (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function make_hamiltonian(J1,D1,N,Nhilbert,irght,ilft,ihfbit,list_1,list_ja,list_jb)\n",
    "    listki = ones(Int,(N+1)*Nhilbert)\n",
    "    loc = ones(Int,(N+1)*Nhilbert)\n",
    "    elemnt = zeros((N+1)*Nhilbert)\n",
    "    for k in 0:N\n",
    "        for i in 0:Nhilbert-1\n",
    "            listki[k*Nhilbert+i+1] = i+1\n",
    "        end\n",
    "    end\n",
    "    for k in 0:N-1\n",
    "        isite1 = k\n",
    "        isite2 = mod((k+1),N)\n",
    "        is1 = 1<<isite1\n",
    "        is2 = 1<<isite2\n",
    "        is0 = is1 + is2\n",
    "        wght = -2.0*J1[k+1]\n",
    "        diag = wght*0.5*D1[k+1]\n",
    "        for i in 0:Nhilbert-1\n",
    "            ii = list_1[i+1]\n",
    "            ibit = ii & is0\n",
    "            if (ibit==0 || ibit==is0)\n",
    "                elemnt[N*Nhilbert+i+1] -= diag\n",
    "                loc[N*Nhilbert+i+1] = i+1\n",
    "            else\n",
    "                elemnt[N*Nhilbert+i+1] += diag\n",
    "                loc[N*Nhilbert+i+1] = i+1\n",
    "                iexchg = xor(ii,is0)\n",
    "                newcfg = get_ja_plus_jb(iexchg,irght,ilft,ihfbit,list_ja,list_jb)\n",
    "                elemnt[k*Nhilbert+i+1] = -wght\n",
    "                loc[k*Nhilbert+i+1] = newcfg+1\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    HamCSC = dropzeros(sparse(listki,loc,elemnt))\n",
    "    return HamCSC\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N=16\n",
      "Sz=0\n",
      "Nup=8\n",
      "Nhilbert=12870\n",
      "ihfbit=256\n",
      "irght,binirght=255 0000000011111111\n",
      "ilft,binilft=65280 1111111100000000\n",
      "iup,biniup=255 0000000011111111\n",
      "  0.045965 seconds (39.98 k allocations: 2.176 MiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([255, 383, 447, 479, 495, 503, 507, 509, 510, 639  …  64896, 65025, 65026, 65028, 65032, 65040, 65056, 65088, 65152, 65280], [0, 0, 1, 0, 2, 1, 2, 0, 3, 3  …  24, 4, 55, 25, 26, 5, 27, 6, 7, 0], [0, 1, 9, 17, 45, 53, 81, 109, 165, 173  …  12669, 12697, 12705, 12761, 12789, 12817, 12825, 12853, 12861, 12869])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "using Arpack\n",
    "using SparseArrays\n",
    "\n",
    "N = 16 # should be N>=4\n",
    "Sz = 0\n",
    "\n",
    "Nup, Nhilbert, ihfbit, irght, ilft, iup = init_parameters(N,Sz)\n",
    "\n",
    "binirght = string(irght,base=2,pad=N)\n",
    "binilft = string(ilft,base=2,pad=N)\n",
    "biniup = string(iup,base=2,pad=N)\n",
    "\n",
    "println(\"N=\",N)\n",
    "println(\"Sz=\",Sz)\n",
    "println(\"Nup=\",Nup)\n",
    "println(\"Nhilbert=\",Nhilbert)\n",
    "println(\"ihfbit=\",ihfbit)\n",
    "println(\"irght,binirght=\",irght,\" \",binirght)\n",
    "println(\"ilft,binilft=\",ilft,\" \",binilft)\n",
    "println(\"iup,biniup=\",iup,\" \",biniup)\n",
    "\n",
    "@time list_1, list_ja, list_jb = make_list(N,Nup,Nhilbert,ihfbit,irght,ilft,iup)\n",
    "\n",
    "#println(\"list_1=\",list_1)\n",
    "#println(\"list_ja=\",list_ja)\n",
    "#println(\"list_jb=\",list_jb)\n",
    "\n",
    "#println(\"\")\n",
    "#println(\"i ii binii ja+jb\")\n",
    "#for i in 0:Nhilbert-1\n",
    "#    ii = list_1[i+1]\n",
    "#    binii = string(ii,base=2,pad=N)\n",
    "#    ind = get_ja_plus_jb(ii,irght,ilft,ihfbit,list_ja,list_jb)\n",
    "#    println(i,\" \",ii,\" \",binii,\" \",ind)\n",
    "#end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.212200 seconds (463.30 k allocations: 35.633 MiB, 2.44% gc time)\n",
      "  2.454857 seconds (7.13 M allocations: 350.538 MiB, 8.14% gc time)\n",
      "# energy:-28.569185442467138 -27.48842671346589 -26.78618970637527 -26.093628229524846 -26.09362822952474\n"
     ]
    }
   ],
   "source": [
    "J1 = ones(N) # J_{ij}>0: AF\n",
    "D1 = ones(N) # D_{ij}>0: AF\n",
    "\n",
    "@time Ham = make_hamiltonian(J1,D1,N,Nhilbert,irght,ilft,ihfbit,list_1,list_ja,list_jb)\n",
    "#println(Ham)\n",
    "\n",
    "@time ene,vec = eigs(Ham, nev=5, which=:LM)\n",
    "#println(\"# GS energy:\",ene[1])\n",
    "println(\"# energy:\",ene[1],\" \",ene[2],\" \",ene[3],\" \",ene[4],\" \",ene[5])\n",
    "#vec_sgn = sign(argmax(vec[:,1]))\n",
    "#println(\"# GS wave function:\")\n",
    "#for i in 0:Nhilbert-1\n",
    "#    ii = list_1[i+1]\n",
    "#    binii = string(ii,base=2,pad=N)\n",
    "#    println(i,\" \",vec[i+1,1]*vec_sgn,\" \",binii)\n",
    "#end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.4",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}