{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lattice Solver Tutorial\n",
    "\n",
    "This tutorial using triqs lattice tools and the Hartree-Fock LatticeSolver for the case of a 2 degenerate bands with local Kanamori interactions at fixed density. At the end, we solve the same problem using the ImpuritySolver with the DMFT loop and arrive at the same solution.\n",
    "\n",
    "We start by importing the required modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: could not identify MPI environment!\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Starting serial run at: 2026-06-12 18:13:06.623061\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from triqs_hartree_fock import ImpuritySolver, LatticeSolver\n",
    "from triqs_hartree_fock.utils import flatten\n",
    "from triqs.gfs import *\n",
    "from triqs.mesh import MeshImFreq\n",
    "from triqs.operators import *\n",
    "from triqs.operators.util import *\n",
    "from h5 import HDFArchive\n",
    "from triqs.lattice.tight_binding import TBLattice\n",
    "from triqs.sumk import *\n",
    "from triqs.lattice import *\n",
    "from scipy.optimize import brentq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we define a 2D tight binding lattice with nearest and next-nearest neighbor hopping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.array([[1, 0.2], [0.2, 1]])\n",
    "tp = np.array([[0.1, 0.05], [0.05, 0.1]])\n",
    "hoppings= {  (1,0)  : t,       \n",
    "        (-1,0) :  t,     \n",
    "        (0,1)  :  t,\n",
    "        (0,-1) :  t,\n",
    "        (1,1)  :  tp,\n",
    "        (-1,-1):  tp,\n",
    "        (1,-1) :  tp,\n",
    "        (-1,1) :  tp}\n",
    "\n",
    "TBL = TBLattice(units=[(1, 0, 0) , (0, 1, 0)], hoppings=hoppings, orbital_positions=[(0., 0., 0.)]*2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define the single particle dispersion h0_k and the other parameters to pass to the solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nk = 10\n",
    "beta = 40\n",
    "U = 3\n",
    "J = 0.5\n",
    "gf_struct = [('up', 2), ('down', 2)]\n",
    "N_target = 2\n",
    "\n",
    "#create a kanamori interaction Hamiltonian\n",
    "Umat, Upmat = U_matrix_kanamori(n_orb=2, U_int=U, J_hund=J)\n",
    "h_int = h_int_kanamori(spin_names=['up', 'down'], n_orb=2, U=Umat, Uprime=Upmat, J_hund=J, off_diag=True)\n",
    "\n",
    "BL = BravaisLattice(units = [(1,0,0) , (0,1,0)])\n",
    "BZ = BrillouinZone(BL)\n",
    "mk = MeshBrZone(BZ, nk)\n",
    "ekup = Gf(mesh=mk, target_shape=[2,2])\n",
    "ekdn = Gf(mesh=mk, target_shape=[2,2])\n",
    "ekup << TBL.fourier(mk)\n",
    "ekdn << TBL.fourier(mk) \n",
    "h0_k = BlockGf(name_list=['up', 'down'], block_list=(ekup, ekdn)) #single particle dispersion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to initialize the solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = LatticeSolver(h0_k=h0_k, gf_struct=gf_struct, beta=beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and finally we call the solve method at fixed density of 2 electrons per site:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "╔╦╗╦═╗╦╔═╗ ╔═╗  ┬ ┬┌─┐\n",
      " ║ ╠╦╝║║═╬╗╚═╗  ├─┤├┤\n",
      " ╩ ╩╚═╩╚═╝╚╚═╝  ┴ ┴└\n",
      "TRIQS: Hartree-Fock solver\n",
      "\n",
      "Running Lattice Solver at fixed density of 2.0000\n",
      "beta = 40.0000\n",
      "h_int =\n",
      "1.5*c_dag('down',0)*c_dag('down',1)*c('down',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',0)*c('up',1)*c('down',1) + 3*c_dag('down',0)*c_dag('up',0)*c('up',0)*c('down',0) + 2*c_dag('down',0)*c_dag('up',1)*c('up',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',1)*c('up',0)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',0)*c('up',1)*c('down',0) + 2*c_dag('down',1)*c_dag('up',0)*c('up',0)*c('down',1) + 3*c_dag('down',1)*c_dag('up',1)*c('up',1)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',1)*c('up',0)*c('down',0) + 1.5*c_dag('up',0)*c_dag('up',1)*c('up',1)*c('up',0)\n",
      "mode: self-consistent\n",
      "Including Fock terms:\n",
      "True\n",
      "Self Consistent Hartree-Fock converged successfully\n",
      "Sigma_HF['up']:\n",
      "[[ 3.25  +0.j -0.0122+0.j]\n",
      " [-0.0122+0.j  3.25  +0.j]]\n",
      "Sigma_HF['down']:\n",
      "[[ 3.25  +0.j -0.0122+0.j]\n",
      " [-0.0122+0.j  3.25  +0.j]]\n",
      "mu = 3.1380\n"
     ]
    }
   ],
   "source": [
    "S.solve(h_int=h_int, N_target=N_target)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DMFT Loop with Impurity Solver:\n",
    "\n",
    "We can arrive at the same solution of the lattice solver using a DMFT loop and impurity solver. However, it's slower than using the lattice solver. We pass a symmetry function that forces the two spins to be equal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu =  -0.11655108826789595\n",
      "\n",
      "╔╦╗╦═╗╦╔═╗ ╔═╗  ┬ ┬┌─┐\n",
      " ║ ╠╦╝║║═╬╗╚═╗  ├─┤├┤\n",
      " ╩ ╩╚═╩╚═╝╚╚═╝  ┴ ┴└\n",
      "TRIQS: Hartree-Fock solver\n",
      "\n",
      "Running Impurity Solver\n",
      "beta = 40.0000\n",
      "h_int =\n",
      "1.5*c_dag('down',0)*c_dag('down',1)*c('down',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',0)*c('up',1)*c('down',1) + 3*c_dag('down',0)*c_dag('up',0)*c('up',0)*c('down',0) + 2*c_dag('down',0)*c_dag('up',1)*c('up',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',1)*c('up',0)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',0)*c('up',1)*c('down',0) + 2*c_dag('down',1)*c_dag('up',0)*c('up',0)*c('down',1) + 3*c_dag('down',1)*c_dag('up',1)*c('up',1)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',1)*c('up',0)*c('down',0) + 1.5*c_dag('up',0)*c_dag('up',1)*c('up',1)*c('up',0)\n",
      "HARTREE SOLVER: starting solver in one-shot mode\n",
      "Including Fock terms:\n",
      "True\n",
      "HARTREE SOLVER: Sigma_HF before iterating['up']:\n",
      "[[0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j]]\n",
      "HARTREE SOLVER: Sigma_HF before iterating['down']:\n",
      "[[0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j]]\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "HARTREE SOLVER: Sigma_HF['up']:\n",
      "[[ 3.25 -0.j -0.009-0.j]\n",
      " [-0.009+0.j  3.25 -0.j]]\n",
      "HARTREE SOLVER: Sigma_HF['down']:\n",
      "[[ 3.25 -0.j -0.009-0.j]\n",
      " [-0.009+0.j  3.25 -0.j]]\n",
      "Final G_dens['up']:\n",
      "[[0.0801-0.j 0.0225-0.j]\n",
      " [0.0225+0.j 0.0801+0.j]]\n",
      "Final G_dens['down']:\n",
      "[[0.0801-0.j 0.0225-0.j]\n",
      " [0.0225+0.j 0.0801+0.j]]\n",
      "mu =  3.1370299374957247\n",
      "\n",
      "╔╦╗╦═╗╦╔═╗ ╔═╗  ┬ ┬┌─┐\n",
      " ║ ╠╦╝║║═╬╗╚═╗  ├─┤├┤\n",
      " ╩ ╩╚═╩╚═╝╚╚═╝  ┴ ┴└\n",
      "TRIQS: Hartree-Fock solver\n",
      "\n",
      "Running Impurity Solver\n",
      "beta = 40.0000\n",
      "h_int =\n",
      "1.5*c_dag('down',0)*c_dag('down',1)*c('down',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',0)*c('up',1)*c('down',1) + 3*c_dag('down',0)*c_dag('up',0)*c('up',0)*c('down',0) + 2*c_dag('down',0)*c_dag('up',1)*c('up',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',1)*c('up',0)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',0)*c('up',1)*c('down',0) + 2*c_dag('down',1)*c_dag('up',0)*c('up',0)*c('down',1) + 3*c_dag('down',1)*c_dag('up',1)*c('up',1)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',1)*c('up',0)*c('down',0) + 1.5*c_dag('up',0)*c_dag('up',1)*c('up',1)*c('up',0)\n",
      "HARTREE SOLVER: starting solver in one-shot mode\n",
      "Including Fock terms:\n",
      "True\n",
      "HARTREE SOLVER: Sigma_HF before iterating['up']:\n",
      "[[ 3.25 -0.j -0.009-0.j]\n",
      " [-0.009+0.j  3.25 -0.j]]\n",
      "HARTREE SOLVER: Sigma_HF before iterating['down']:\n",
      "[[ 3.25 -0.j -0.009-0.j]\n",
      " [-0.009+0.j  3.25 -0.j]]\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "HARTREE SOLVER: Sigma_HF['up']:\n",
      "[[ 3.25  +0.j -0.0112-0.j]\n",
      " [-0.0112+0.j  3.25  +0.j]]\n",
      "HARTREE SOLVER: Sigma_HF['down']:\n",
      "[[ 3.25  +0.j -0.0112-0.j]\n",
      " [-0.0112+0.j  3.25  +0.j]]\n",
      "Final G_dens['up']:\n",
      "[[0.4998-0.j 0.0232+0.j]\n",
      " [0.0232-0.j 0.4998-0.j]]\n",
      "Final G_dens['down']:\n",
      "[[0.4998-0.j 0.0232+0.j]\n",
      " [0.0232-0.j 0.4998-0.j]]\n",
      "mu =  3.1377229169767435\n",
      "\n",
      "╔╦╗╦═╗╦╔═╗ ╔═╗  ┬ ┬┌─┐\n",
      " ║ ╠╦╝║║═╬╗╚═╗  ├─┤├┤\n",
      " ╩ ╩╚═╩╚═╝╚╚═╝  ┴ ┴└\n",
      "TRIQS: Hartree-Fock solver\n",
      "\n",
      "Running Impurity Solver\n",
      "beta = 40.0000\n",
      "h_int =\n",
      "1.5*c_dag('down',0)*c_dag('down',1)*c('down',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',0)*c('up',1)*c('down',1) + 3*c_dag('down',0)*c_dag('up',0)*c('up',0)*c('down',0) + 2*c_dag('down',0)*c_dag('up',1)*c('up',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',1)*c('up',0)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',0)*c('up',1)*c('down',0) + 2*c_dag('down',1)*c_dag('up',0)*c('up',0)*c('down',1) + 3*c_dag('down',1)*c_dag('up',1)*c('up',1)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',1)*c('up',0)*c('down',0) + 1.5*c_dag('up',0)*c_dag('up',1)*c('up',1)*c('up',0)\n",
      "HARTREE SOLVER: starting solver in one-shot mode\n",
      "Including Fock terms:\n",
      "True\n",
      "HARTREE SOLVER: Sigma_HF before iterating['up']:\n",
      "[[ 3.25  +0.j -0.0112-0.j]\n",
      " [-0.0112+0.j  3.25  +0.j]]\n",
      "HARTREE SOLVER: Sigma_HF before iterating['down']:\n",
      "[[ 3.25  +0.j -0.0112-0.j]\n",
      " [-0.0112+0.j  3.25  +0.j]]\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "HARTREE SOLVER: Sigma_HF['up']:\n",
      "[[ 3.25  -0.j -0.0119+0.j]\n",
      " [-0.0119-0.j  3.25  -0.j]]\n",
      "HARTREE SOLVER: Sigma_HF['down']:\n",
      "[[ 3.25  -0.j -0.0119+0.j]\n",
      " [-0.0119-0.j  3.25  -0.j]]\n",
      "Final G_dens['up']:\n",
      "[[0.4999+0.j 0.024 -0.j]\n",
      " [0.024 +0.j 0.4999+0.j]]\n",
      "Final G_dens['down']:\n",
      "[[0.4999+0.j 0.024 -0.j]\n",
      " [0.024 +0.j 0.4999+0.j]]\n",
      "mu =  3.1379142043949604\n",
      "\n",
      "╔╦╗╦═╗╦╔═╗ ╔═╗  ┬ ┬┌─┐\n",
      " ║ ╠╦╝║║═╬╗╚═╗  ├─┤├┤\n",
      " ╩ ╩╚═╩╚═╝╚╚═╝  ┴ ┴└\n",
      "TRIQS: Hartree-Fock solver\n",
      "\n",
      "Running Impurity Solver\n",
      "beta = 40.0000\n",
      "h_int =\n",
      "1.5*c_dag('down',0)*c_dag('down',1)*c('down',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',0)*c('up',1)*c('down',1) + 3*c_dag('down',0)*c_dag('up',0)*c('up',0)*c('down',0) + 2*c_dag('down',0)*c_dag('up',1)*c('up',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',1)*c('up',0)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',0)*c('up',1)*c('down',0) + 2*c_dag('down',1)*c_dag('up',0)*c('up',0)*c('down',1) + 3*c_dag('down',1)*c_dag('up',1)*c('up',1)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',1)*c('up',0)*c('down',0) + 1.5*c_dag('up',0)*c_dag('up',1)*c('up',1)*c('up',0)\n",
      "HARTREE SOLVER: starting solver in one-shot mode\n",
      "Including Fock terms:\n",
      "True\n",
      "HARTREE SOLVER: Sigma_HF before iterating['up']:\n",
      "[[ 3.25  -0.j -0.0119+0.j]\n",
      " [-0.0119-0.j  3.25  -0.j]]\n",
      "HARTREE SOLVER: Sigma_HF before iterating['down']:\n",
      "[[ 3.25  -0.j -0.0119+0.j]\n",
      " [-0.0119-0.j  3.25  -0.j]]\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "HARTREE SOLVER: Sigma_HF['up']:\n",
      "[[ 3.25  -0.j -0.0121-0.j]\n",
      " [-0.0121+0.j  3.25  -0.j]]\n",
      "HARTREE SOLVER: Sigma_HF['down']:\n",
      "[[ 3.25  -0.j -0.0121-0.j]\n",
      " [-0.0121+0.j  3.25  -0.j]]\n",
      "Final G_dens['up']:\n",
      "[[0.5   -0.j 0.0243-0.j]\n",
      " [0.0243+0.j 0.5   +0.j]]\n",
      "Final G_dens['down']:\n",
      "[[0.5   -0.j 0.0243-0.j]\n",
      " [0.0243+0.j 0.5   +0.j]]\n",
      "mu =  3.1379721665000178\n",
      "\n",
      "╔╦╗╦═╗╦╔═╗ ╔═╗  ┬ ┬┌─┐\n",
      " ║ ╠╦╝║║═╬╗╚═╗  ├─┤├┤\n",
      " ╩ ╩╚═╩╚═╝╚╚═╝  ┴ ┴└\n",
      "TRIQS: Hartree-Fock solver\n",
      "\n",
      "Running Impurity Solver\n",
      "beta = 40.0000\n",
      "h_int =\n",
      "1.5*c_dag('down',0)*c_dag('down',1)*c('down',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',0)*c('up',1)*c('down',1) + 3*c_dag('down',0)*c_dag('up',0)*c('up',0)*c('down',0) + 2*c_dag('down',0)*c_dag('up',1)*c('up',1)*c('down',0) + 0.5*c_dag('down',0)*c_dag('up',1)*c('up',0)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',0)*c('up',1)*c('down',0) + 2*c_dag('down',1)*c_dag('up',0)*c('up',0)*c('down',1) + 3*c_dag('down',1)*c_dag('up',1)*c('up',1)*c('down',1) + 0.5*c_dag('down',1)*c_dag('up',1)*c('up',0)*c('down',0) + 1.5*c_dag('up',0)*c_dag('up',1)*c('up',1)*c('up',0)\n",
      "HARTREE SOLVER: starting solver in one-shot mode\n",
      "Including Fock terms:\n",
      "True\n",
      "HARTREE SOLVER: Sigma_HF before iterating['up']:\n",
      "[[ 3.25  -0.j -0.0121-0.j]\n",
      " [-0.0121+0.j  3.25  -0.j]]\n",
      "HARTREE SOLVER: Sigma_HF before iterating['down']:\n",
      "[[ 3.25  -0.j -0.0121-0.j]\n",
      " [-0.0121+0.j  3.25  -0.j]]\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "\n",
      "HARTREE SOLVER: fixing dc to 0.0000\n",
      "HARTREE SOLVER: Sigma_HF['up']:\n",
      "[[ 3.25  +0.j -0.0122-0.j]\n",
      " [-0.0122+0.j  3.25  +0.j]]\n",
      "HARTREE SOLVER: Sigma_HF['down']:\n",
      "[[ 3.25  +0.j -0.0122-0.j]\n",
      " [-0.0122+0.j  3.25  +0.j]]\n",
      "Final G_dens['up']:\n",
      "[[0.5   -0.j 0.0244+0.j]\n",
      " [0.0244-0.j 0.5   -0.j]]\n",
      "Final G_dens['down']:\n",
      "[[0.5   -0.j 0.0244+0.j]\n",
      " [0.0244-0.j 0.5   -0.j]]\n"
     ]
    }
   ],
   "source": [
    "SK = SumkDiscreteFromLattice(lattice=TBL, n_points=nk)\n",
    "iw_mesh = MeshImFreq(beta=beta, statistic='Fermion', n_iw=500)\n",
    "sigma = Gf(mesh=iw_mesh, target_shape=[2,2])\n",
    "Sigma = BlockGf(name_list=['up', 'down'], block_list=(sigma,sigma), make_copies=True)\n",
    "Gloc = Sigma.copy()\n",
    "N_target = 2\n",
    "mu = 0\n",
    "\n",
    "#symmetry function forcing two spins to be equal\n",
    "def make_spins_equal(Sigma):\n",
    "    Symmetrized_Sigma = Sigma.copy()\n",
    "    Symmetrized_Sigma['up'] = 0.5*(Sigma['up']  + Sigma['down'])\n",
    "    Symmetrized_Sigma['down'] = Symmetrized_Sigma['up']\n",
    "    return Symmetrized_Sigma\n",
    "\n",
    "S_imp = ImpuritySolver(gf_struct=gf_struct, mesh=iw_mesh, symmetries=[make_spins_equal])\n",
    "\n",
    "converged = False\n",
    "while not converged:\n",
    "    for name, bl in gf_struct:\n",
    "        Sigma[name] << S_imp.Sigma_HF[name]\n",
    "    mu = brentq(lambda mu: SK(mu=mu, Sigma=Sigma).total_density().real - N_target, -5, 5)\n",
    "    print('mu = ', mu)\n",
    "    Gloc << SK(mu=mu, Sigma=Sigma)\n",
    "    S_imp.G0_iw << inverse(inverse(Gloc) + Sigma)\n",
    "    Sigma_old = S_imp.Sigma_HF.copy()\n",
    "    S_imp.solve(h_int=h_int)\n",
    "    if np.allclose(flatten(Sigma_old), flatten(S_imp.Sigma_HF), rtol=0, atol=1e-4):\n",
    "        converged = True"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}