Multi-orbital Hartree-Fock calculations in Julia

hf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "upper-dover",
   "metadata": {},
   "source": [
    "# Hartree-Fock calculation of multi-orbital systems\n",
    "\n",
    "We consider the Hamiltonian\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\sum_{ij=1}^N t_{ij} \\hat{c}^\\dagger_i \\hat{c}_j +  \\frac{1}{4} \\sum_{ijkl=1}^N U_{ijlk}\\hat{c}^\\dagger_i \\hat{c}^\\dagger_j \\hat{c}_k \\hat{c}_l,\n",
    "$$\n",
    "where $i$ and $j$ are composite indices of spin and orbitals, $N$ is the number of spin orbitals.\n",
    "The Coulomb utensor $U$ is antisymmetrized as $U_{ijlk} = -U_{ijkl} = -U_{jilk}$.\n",
    "\n",
    "After simple manipulations, we obtain a mean-field Hamiltonian\n",
    "\n",
    "$$\n",
    "\\hat{H}_\\mathrm{MF} = \\sum_{ij=1}^N T_{ij} \\hat{c}_i^\\dagger \\hat{c}_j -\\frac{1}{2}\\sum_{ijkl=1}^N U_{ijlk}\\langle \\hat{c}^\\dagger_i \\hat{c}_l\\rangle \\langle \\hat{c}^\\dagger_j \\hat{c}_k\\rangle,\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "T_{ij} = t_{ij} + \\sum_{kl=1}^N U_{i k j l} \\langle \\hat{c}^\\dagger_k \\hat{c}_l\\rangle.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "neural-drink",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "solve_hf"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "using OMEinsum\n",
    "\n",
    "\"\"\"\n",
    "Compute hopping matrix of the MF Hamiltonian, tmat\n",
    "\"\"\"\n",
    "function compute_tmat(tmat, utensor, rhomat)\n",
    "    tmat = copy(tmat)\n",
    "    N = size(tmat)[1]\n",
    "    for i in 1:N, j in 1:N\n",
    "        for k in 1:N, l in 1:N\n",
    "            tmat[i,j] += utensor[i,k,j,l] * rhomat[k,l]\n",
    "        end\n",
    "    end\n",
    "    tmat\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "Compute density matrix rhomat\n",
    "\"\"\"\n",
    "function compute_rhomat!(rhomat, evecs, nelec)\n",
    "    rhomat .= 0.0\n",
    "    N = size(rhomat)[1]\n",
    "    for e in 1:nelec\n",
    "        for i in 1:N, j in 1:N\n",
    "           rhomat[i, j] += evecs[i,e] * evecs[j,e]\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "Self-consistent calculations\n",
    "\"\"\"\n",
    "function solve_hf(tmat::Array{Float64,2}, utensor::Array{Float64,4}, rhomat0::Array{Float64,2}, nelec, niter, mixing)\n",
    "    rhomat = copy(rhomat0)\n",
    "    rhomat_out = copy(rhomat)\n",
    "    sum_eigen = 0.0\n",
    "    for iter in 1:niter\n",
    "        Tmat = compute_tmat(tmat, utensor, rhomat)\n",
    "        e = eigen(Hermitian(Tmat))\n",
    "        sum_eigen = sum(e.values[1:nelec])\n",
    "        compute_rhomat!(rhomat_out, e.vectors, nelec)\n",
    "        @. rhomat = (1-mixing) * rhomat + mixing * rhomat_out\n",
    "    end\n",
    "    # Compute energy\n",
    "    ene = sum_eigen - 0.5 * ein\"ijlk,il,jk->\"(utensor, rhomat, rhomat)[1]\n",
    "    rhomat, ene\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "isolated-closer",
   "metadata": {},
   "source": [
    "We consider the two-site Hubbard model with onsite repulsion $U$.\n",
    "For each site, the Coulomb interaction can be represented as\n",
    "\n",
    "$$\n",
    "U n_{\\uparrow} n_{\\downarrow} = \\frac{1}{4} \\sum_{s_1 s_2 s_3 s_4} U_{s_1 s_2 s_3 s_4} \\hat{c}^\\dagger_{s_1} \\hat{c}^\\dagger_{s_2} c_{s_4} c_{s_3},\n",
    "$$\n",
    "\n",
    "where $U_{\\uparrow\\downarrow\\uparrow\\downarrow}= U_{\\downarrow\\uparrow\\downarrow\\uparrow}=U$ and $U_{\\uparrow\\downarrow\\downarrow\\uparrow}= \\bar{U}_{\\downarrow\\uparrow\\uparrow\\downarrow}=-U$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "large-currency",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nsites = 2\n",
    "up = 1\n",
    "dn = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fitting-performer",
   "metadata": {},
   "source": [
    "## Construct Coulomb utensor\n",
    "Each spin-orbital index runs as (site1, up), (site1, dn), (site2, up), (site2, dn)\n",
    "because Julia uses column major for memory layout. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "viral-headline",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "generate_onsiteU (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function generate_onsiteU(nsites, U)\n",
    "    utensor = zeros(Float64, 2, nsites, 2, nsites, 2, nsites, 2, nsites)\n",
    "    for isite in 1:nsites\n",
    "        utensor[up, isite, dn, isite, up, isite, dn, isite] = U\n",
    "        utensor[dn, isite, up, isite, dn, isite, up, isite] = U\n",
    "        utensor[up, isite, dn, isite, dn, isite, up, isite] = -U\n",
    "        utensor[dn, isite, up, isite, up, isite, dn, isite] = -U\n",
    "    end\n",
    "    reshape(utensor, (2*nsites, 2*nsites, 2*nsites, 2*nsites))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "distinct-superintendent",
   "metadata": {},
   "source": [
    "## Construct hopping matrix (spin-diagonal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "underlying-nudist",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "generate_NNhopping (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function generate_NNhopping(nsites, t)\n",
    "    tmat = zeros(Float64, 2, nsites, 2, nsites)\n",
    "    for isite in 1:nsites-1\n",
    "        for spin in 1:2\n",
    "            tmat[spin, isite, spin, isite+1] = t\n",
    "            tmat[spin, isite+1, spin, isite] = t\n",
    "        end\n",
    "    end\n",
    "    reshape(tmat, 2*nsites, 2*nsites)\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accomplished-overview",
   "metadata": {},
   "source": [
    "## Initial density matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "yellow-integrity",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "generate_rhomat (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "generate_rhomat(occup) = diagm(occup)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "buried-upper",
   "metadata": {},
   "source": [
    "## Tests (isolated & fully occupied case)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "static-wildlife",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of electrons = 4.0\n",
      "  6.685071 seconds (24.60 M allocations: 1.219 GiB, 6.74% gc time)\n"
     ]
    }
   ],
   "source": [
    "U = 4.0\n",
    "t = 0.0\n",
    "nelec = 4\n",
    "\n",
    "niter = 100\n",
    "tmat = generate_NNhopping(nsites, 0.0)\n",
    "utensor = generate_onsiteU(nsites, U)\n",
    "rhomat0 = generate_rhomat([1.0, 1.0, 1.0, 1.0])\n",
    "@assert sum(diag(rhomat0)) ≈ 4\n",
    "println(\"Total number of electrons = \", sum(diag(rhomat0)))\n",
    "rhomat, ene = @time solve_hf(tmat, utensor, rhomat0, 4, niter, 0.1)\n",
    "@assert diag(rhomat0) ≈ ones(4)\n",
    "@assert ene ≈ U * nsites"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acceptable-resident",
   "metadata": {},
   "source": [
    "## Self-consistent calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "current-selling",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of electrons = 2.0\n",
      "  0.000770 seconds (1.63 k allocations: 312.422 KiB)\n"
     ]
    }
   ],
   "source": [
    "U = 4.0\n",
    "t = 1.0\n",
    "nelec = 2\n",
    "\n",
    "niter = 100\n",
    "tmat = generate_NNhopping(nsites, t)\n",
    "utensor = generate_onsiteU(nsites, U)\n",
    "rhomat0 = generate_rhomat([1.0, 0.0, 0.0, 1.0])\n",
    "@assert sum(diag(rhomat0)) ≈ nelec\n",
    "println(\"Total number of electrons = \", sum(diag(rhomat0)))\n",
    "rhomat, ene = @time solve_hf(tmat, utensor, rhomat0, nelec, niter, 0.1)\n",
    ";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "southern-bridge",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.500014617187124"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "offensive-prison",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "after-launch",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plastic-personality",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.3",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}