Supplemental_II.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix B: Non-uniqueness of mean-field projection and of the DMET fixed point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.Non-uniqueness of the mean-field projection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We explicitly show non-uniqueness of the approximate (mean-field) projection by constructing two non-interacting system that have the same impurity 1RDM $\\gamma^{s}_{\\rm imp}(i,j)$ but different orbitals and thus projections. \n",
    "\n",
    "As our example system, we choose a translation invariant Hubbard system with periodic boundary conditions, occupied by 3 particles in the ground state. Further, we choose our impurity to consist of the first two lattice sites of this system.<br>\n",
    "Our Hamiltonian thus reads:\n",
    " \\begin{align}\n",
    "    \\boldsymbol{h}^{s}=\\begin{pmatrix}\n",
    "0& -1    &  0 & 0& 0 &-1\\\\\n",
    "-1& 0   &  -1 & 0& 0 &0\\\\\n",
    "0& -1   &  0 & -1 & 0 &0\\\\\n",
    "0& 0   &  -1 & 0 & -1 &0\\\\\n",
    "0& 0   &  0 & -1 & 0 &-1\\\\\n",
    "-1& 0   &  0 & 0 & -1 &0\\\\\n",
    "\\end{pmatrix}\n",
    "\\label{eq:ham_targ_6}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "from numpy.linalg import multi_dot\n",
    "np.set_printoptions(precision=5)\n",
    "import copy\n",
    "\n",
    "\n",
    "h = -np.ones((6,6))\n",
    "h *= np.tri(*h.shape, k=1)\n",
    "np.fill_diagonal(h,0)\n",
    "h *= 1 - np.tri(*h.shape, k=-2)\n",
    "h[0,5]=-1\n",
    "h[5,0]=-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0. -1. -0. -0. -0. -1.]\n",
      " [-1.  0. -1. -0. -0. -0.]\n",
      " [-0. -1.  0. -1. -0. -0.]\n",
      " [-0. -0. -1.  0. -1. -0.]\n",
      " [-0. -0. -0. -1.  0. -1.]\n",
      " [-1. -0. -0. -0. -1.  0.]]\n"
     ]
    }
   ],
   "source": [
    "print(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Diagonalizing $\\boldsymbol{h}^{s}$ leads to six eigenstates $\\{\\phi_{1}, ..., \\phi_{6}\\}$. <br>\n",
    "We occupy the three lowest eigenstates $\\{\\phi_{1},\\phi_{2},\\phi_{3}\\}$ with ground-state energy $E=-4$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "epsilon, phi = LA.eigh(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-4.0\n"
     ]
    }
   ],
   "source": [
    "E=epsilon[0]+epsilon[1]+epsilon[2]\n",
    "print(np.around(E,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 1RDM of  $\\{\\phi_{1},\\phi_{2}, \\phi_{3}\\}$ reads:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \\begin{align}\n",
    "    \\boldsymbol{\\gamma}^{s}=\\begin{pmatrix}\n",
    "0.5& 0.33333    &  0 & -0.16666& 0 &0.33333\\\\\n",
    "0.33333& 0.5   &  0.33333 & 0& -0.16666 &0\\\\\n",
    "0& 0.33333   &  0.5 & 0.33333 & 0 &-0.16666\\\\\n",
    "-0.16666& 0   &  0.33333 & 0.5 & 0.33333 &0\\\\\n",
    "0& -0.166666  &  0 & 0.33333& 0.5 &0.333333\\\\\n",
    "0.333333& 0   &  -0.166666 & 0 & 0.333333 &0.5\\\\\n",
    "\\end{pmatrix}\n",
    "\\end{align}."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 5.00000e-01  3.33333e-01  1.60448e-16 -1.66667e-01  1.47285e-16\n",
      "   3.33333e-01]\n",
      " [ 3.33333e-01  5.00000e-01  3.33333e-01 -2.30550e-17 -1.66667e-01\n",
      "   8.26456e-17]\n",
      " [ 1.60448e-16  3.33333e-01  5.00000e-01  3.33333e-01  6.40920e-17\n",
      "  -1.66667e-01]\n",
      " [-1.66667e-01 -2.30550e-17  3.33333e-01  5.00000e-01  3.33333e-01\n",
      "  -2.12794e-17]\n",
      " [ 1.47285e-16 -1.66667e-01  6.40920e-17  3.33333e-01  5.00000e-01\n",
      "   3.33333e-01]\n",
      " [ 3.33333e-01  8.26456e-17 -1.66667e-01 -2.12794e-17  3.33333e-01\n",
      "   5.00000e-01]]\n"
     ]
    }
   ],
   "source": [
    "gamma_s=np.dot(phi[:,0:3],phi[:,0:3].transpose())\n",
    "print(gamma_s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting impurity (which are the first two lattice sites) 1RDM on $A \\equiv\\{1,2 \\}$ is then\n",
    "\\begin{align}\n",
    "    \\boldsymbol{\\gamma}^{s}_{\\rm imp}=\\begin{pmatrix}\n",
    "0.5& 0.33333    \\\\\n",
    "0.33333& 0.5    \\\\\n",
    "\\end{pmatrix}\n",
    "\\label{eq:gamma_t_imp}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.5     0.33333]\n",
      " [0.33333 0.5    ]]\n"
     ]
    }
   ],
   "source": [
    "gamma_imp=gamma_s[0:2,0:2]\n",
    "print(gamma_imp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to construct a different system with the same impurity 1RDM as its three-particle ground-state we first \n",
    "diagonalize $\\boldsymbol{\\gamma}_{\\rm imp}$ and get the eigenvalues and the eigenvectors of the impurity 1RDM as\n",
    "\\begin{align}\n",
    "n_1^{\\rm imp}&=0.16666 \\equiv \\| \\tilde{\\varphi}_1\\|_A^2\\nonumber\\\\\n",
    "\\varphi^{A}_1 &=  0.70711 |1\\rangle\n",
    " -0.70711  |2\\rangle.\n",
    "\\end{align}\n",
    "\\begin{align}\n",
    "n_2^{\\rm imp}&=0.83333 \\equiv \\|\\tilde{\\varphi}_2\\|_A^2 \\nonumber\\\\\n",
    "\\varphi^{A}_2&= 0.70711 |1\\rangle\n",
    " +0.70711  |2\\rangle\n",
    "  \\label{eq:imp_occ_2}\n",
    "\\end{align}  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1666666666666666\n",
      "[ 0.70711 -0.70711]\n",
      "0.833333333333333\n",
      "[-0.70711 -0.70711]\n"
     ]
    }
   ],
   "source": [
    "n_imp, phi_imp=LA.eigh(gamma_imp)\n",
    "for i in range (0,2):\n",
    "    print(n_imp[i])\n",
    "    print(phi_imp[:,i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the whole system consists of 6 lattice sites and we choose our impurity to consist of 2 lattice sites, the environment $B$ is four-dimensional, thus consisting of four basis functions.<br> \n",
    "We can choose problem adopted ones by just diagonalizing the environment 1RDM of the original $\\gamma^{s}(i,j)$ and use two of them to build our CAS space, i.e."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 5.00000e-01  3.33333e-01  6.40920e-17 -1.66667e-01]\n",
      " [ 3.33333e-01  5.00000e-01  3.33333e-01 -2.12794e-17]\n",
      " [ 6.40920e-17  3.33333e-01  5.00000e-01  3.33333e-01]\n",
      " [-1.66667e-01 -2.12794e-17  3.33333e-01  5.00000e-01]]\n"
     ]
    }
   ],
   "source": [
    "gamma_env=gamma_s[2:6,2:6]\n",
    "print(gamma_env)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This leads to\n",
    "  \\begin{align}\n",
    "n_1&=0\\equiv \\|\\tilde{\\varphi}_1\\|_B^2\\\\\n",
    "\\varphi^{B}_1&= - 0.31623|3\\rangle+0.63246|4\\rangle -0.63246|5\\rangle +0.31623|6\\rangle \n",
    "  \\label{eq:env_unocc}\n",
    "\\end{align}\n",
    "\n",
    "  \\begin{align}\n",
    "n_2&=  0.16666 \\equiv \\|\\tilde{\\varphi}_2\\|_B^2\\\\   \n",
    "\\varphi^{B}_2&=\n",
    "  0.63246 |3\\rangle     \n",
    " -0.31623 |4\\rangle\n",
    " -0.31623 |5 \\rangle   \n",
    "  +0.63246 |6\\rangle \\nonumber\n",
    "  \\label{eq:n_2_env}\n",
    "  \\end{align}\n",
    "   \\begin{align}\n",
    "  n_3&= 0.83333 \\equiv \\|\\tilde{\\varphi}_3\\|_B^2 \\\\\n",
    "  \\varphi^{B}_3&=\n",
    " -0.63246|3\\rangle    \n",
    " -0.31623 |4\\rangle+\n",
    "  0.31623 |5\\rangle+     \n",
    "  0.63246 |6\\rangle   \n",
    "  \\label{eq:n_3_env}\n",
    "  \\end{align}\n",
    " \\begin{align}\n",
    " \\label{eq:phi_core}\n",
    " n_4&= 1.0 \\equiv \\| \\tilde{\\varphi}_4\\|_B^2\\nonumber\\\\\n",
    " \\varphi^{B}_4&=\n",
    "  0.31623|3\\rangle     + 0.63246  |4\\rangle\n",
    "+ 0.63246   |5\\rangle+ 0.31623 |6\\rangle.  \n",
    "  \\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.8219613577756638e-16\n",
      "[-0.31623  0.63246 -0.63246  0.31623]\n",
      "0.16666666666666677\n",
      "[ 0.63246 -0.31623 -0.31623  0.63246]\n",
      "0.8333333333333327\n",
      "[ 0.63246  0.31623 -0.31623 -0.63246]\n",
      "0.9999999999999998\n",
      "[-0.31623 -0.63246 -0.63246 -0.31623]\n"
     ]
    }
   ],
   "source": [
    "n_env, phi_env = LA.eigh(gamma_env)\n",
    "for i in range (0,4):\n",
    "    print(n_env[i])\n",
    "    print(phi_env[:,i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can construct the CAS orbitals that would be used in the auxiliary projection of the target Hamiltonian (where for the purpose of the example is not interacting but in a real application one would be interested in interacting Hamiltonians). The first two CAS orbitals can be always chosen as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "    \\varphi^{\\rm CAS}_1&= |1\\rangle, \n",
    "  \\label{eq:phi_1_cas}\n",
    "  \\\\\n",
    "    \\varphi^{\\rm CAS}_2&= |2\\rangle.\n",
    "  \\label{eq:phi_2_cas}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we could have used the two orbitals $\\phi^A$ of the impurity submatrix as we have discussed in Appendix A.<br>\n",
    "Because the fourth eigenvector of the environment submatrix is discarded in the usual approximate projection (unentangled occupied/core orbital) and the first orbital is perpendicular to the subspace of the three lowest orbitals, we build the other CAS (environmental) orbitals from the remaining orbitals as\n",
    "\\begin{align}\n",
    "    \\varphi^{\\rm CAS}_3&=\\varphi^{B}_2,\n",
    "  \\label{eq:phi_3_cas}\n",
    "\\\\\n",
    "  \\varphi^{\\rm CAS}_4&=\\varphi^{B}_3.\n",
    "  \\label{eq:phi_4_cas}\n",
    "\\end{align}\n",
    "While we do not need this CAS orbitals in this section, they will become important in the next. Further, they will show that we get a very different projection when we compare to the CAS from the different Hamiltonian that we construct next.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.       0.       0.       0.     ]\n",
      " [ 0.       1.       0.       0.     ]\n",
      " [ 0.       0.       0.63246  0.63246]\n",
      " [ 0.       0.      -0.31623  0.31623]\n",
      " [ 0.       0.      -0.31623 -0.31623]\n",
      " [ 0.       0.       0.63246 -0.63246]]\n"
     ]
    }
   ],
   "source": [
    "phi_CAS=np.zeros((6,4))\n",
    "for i in range(2):\n",
    "    phi_CAS[i,i]=1\n",
    "    phi_CAS[2:6,i+2]=phi_env[:,i+1]\n",
    "print(phi_CAS)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we now take $\\varphi_3^{B}$ and $\\varphi_4^{B}$ and define\n",
    "\\begin{align}\n",
    "    \\tilde{\\varphi}'_1 &= \\sqrt{n_1^{\\rm imp}} \\varphi^{A}_1 + \\underbrace{\\sqrt{1-n_1^{\\rm imp}}}_{=\\sqrt{0.83333}} \\varphi_1^B\n",
    "    \\\\\n",
    "    \\tilde{\\varphi}'_2 &= \\sqrt{n_2^{\\rm imp}} \\varphi^{A}_2 + \\underbrace{\\sqrt{1-n_2^{\\rm imp}}}_{=\\sqrt{0.16666}} \\varphi_4^B\n",
    "\\end{align}\n",
    "as well as\n",
    "\\begin{align}\n",
    "    \\tilde{\\varphi}'_3 &= \\varphi^B_2\\\\\n",
    "    \\tilde{\\varphi}'_4 &= \\varphi^B_1.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.28868 -0.6455   0.       0.     ]\n",
      " [-0.28868 -0.6455   0.       0.     ]\n",
      " [-0.28868 -0.1291   0.63246  0.63246]\n",
      " [ 0.57735 -0.2582  -0.31623  0.31623]\n",
      " [-0.57735 -0.2582  -0.31623 -0.31623]\n",
      " [ 0.28868 -0.1291   0.63246 -0.63246]]\n"
     ]
    }
   ],
   "source": [
    "phi_tilde=np.zeros((6,4))\n",
    "\n",
    "phi_tilde[0:2,0]=n_imp[0]**0.5*phi_imp[0:2,0]\n",
    "phi_tilde[2:6,0]=(1-n_imp[0])**0.5*phi_env[0:4,0]\n",
    "phi_tilde[0:2,1]=n_imp[1]**0.5*phi_imp[0:2,1]\n",
    "phi_tilde[2:6,1]=(1-n_imp[1])**0.5*phi_env[0:4,3]\n",
    "\n",
    "phi_tilde[2:6,2]=phi_env[:,1]\n",
    "phi_tilde[2:6,3]=phi_env[:,2]\n",
    "print(phi_tilde)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have now four orthogonal vectors we would still need to choose two orthonormal ones to fill up all of the six dimensional space. However, since we only want to construct a Hamiltonian that has the same impurity 1RDM in the ground-state three-particle sector we leave them undefined but instead choose a set of random numbers $\\epsilon_1' \\leq \\epsilon'_2 \\leq \\epsilon'_3 <  \\epsilon'_4 \\leq \\epsilon'_5 = \\epsilon'_6 =0$. For definiteness, let us choose $\\epsilon_1' = -4$, $\\epsilon'_2= -3 $, $\\epsilon'_3=-2$, $\\epsilon'_4= -1$.\n",
    "With these eigenvalues $\\epsilon'_i$ and the $\\tilde{\\varphi}'_i$ we can build a new Hamiltonian   $\\boldsymbol{h}^{s}_{\\rm new}=\\sum_{\\mu =1}^{4} \\epsilon'_{\\mu} \\tilde{\\varphi}'_{\\mu}(i)\\tilde{\\varphi}'^{*}_{\\mu (j)}$:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.58333 -0.91667  0.08333 -1.16667  0.16667 -0.58333]\n",
      " [-0.91667 -1.58333 -0.58333  0.16667 -1.16667  0.08333]\n",
      " [ 0.08333 -0.58333 -1.58333  0.76667 -0.16667 -0.11667]\n",
      " [-1.16667  0.16667  0.76667 -1.83333  1.03333 -0.16667]\n",
      " [ 0.16667 -1.16667 -0.16667  1.03333 -1.83333  0.76667]\n",
      " [-0.58333  0.08333 -0.11667 -0.16667  0.76667 -1.58333]]\n"
     ]
    }
   ],
   "source": [
    "epsilon_prime=np.zeros(4)\n",
    "h_new_s=np.zeros((6,6))\n",
    "for i in range(0,4):\n",
    "    epsilon_prime[i]=-4+i   \n",
    "\n",
    "h_new_s=np.zeros((6,6))\n",
    "h_new_s=multi_dot([phi_tilde,np.diag(epsilon_prime),phi_tilde.transpose()])\n",
    "    \n",
    "print(h_new_s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we diagonalize the Hamiltonian and take the three lowest eigenvectors (which would be identical to $\\tilde{\\varphi}'_i$ for i=1,2,3 since we have build this Hamiltonian to have these as eigenvectors) we find the three-particle ground-state 1RDM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 5.00000e-01  3.33333e-01  2.55739e-16  3.33333e-01  7.44393e-17\n",
      "   1.66667e-01]\n",
      " [ 3.33333e-01  5.00000e-01  1.66667e-01  7.35290e-17  3.33333e-01\n",
      "  -1.45922e-16]\n",
      " [ 2.55739e-16  1.66667e-01  5.00000e-01 -3.33333e-01  3.24340e-16\n",
      "   3.33333e-01]\n",
      " [ 3.33333e-01  7.35290e-17 -3.33333e-01  5.00000e-01 -1.66667e-01\n",
      "   4.75155e-17]\n",
      " [ 7.44393e-17  3.33333e-01  3.24340e-16 -1.66667e-01  5.00000e-01\n",
      "  -3.33333e-01]\n",
      " [ 1.66667e-01 -1.45922e-16  3.33333e-01  4.75155e-17 -3.33333e-01\n",
      "   5.00000e-01]]\n"
     ]
    }
   ],
   "source": [
    "gamma_s_new=np.dot(phi_tilde[:,0:3],phi_tilde[:,0:3].transpose())\n",
    "print(gamma_s_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the corresponding ground-state energy is $E'=\\epsilon'_1+\\epsilon'_2+\\epsilon'_3 = -9$. By construction the 1RDM agrees on the impurity but the rest is different. Also, the CAS orbitals that are used in the projection will be different. In our case they become (besides the first two that are always the same)\n",
    "\\begin{align}\n",
    "    \\varphi'^{\\rm CAS}_3&=\\varphi^{B}_4,\n",
    "  \\label{eq:phi_3_cas_prime}\n",
    "\\\\\n",
    "  \\varphi'^{\\rm CAS}_4&=\\varphi^{B}_1.\n",
    "  \\label{eq:phi_4_cas_prime}\n",
    "\\end{align}\n",
    " However, $\\varphi^{B}_4$ and $\\varphi^{B}_1$ were orbitals that did not belong to the original CAS. That means that also the projection constructed from the hamiltonian $\\boldsymbol{h}^s_{\\rm new}$ will look very different from the one of $\\boldsymbol{h}^{s}$. Thus in this example we have highlighted that by the requirement that the impurity 1RDM is the same there is the possibility to construct completely different projections even in a very simple setting where also the target system is non-interacting and we just consider only a few sites.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Non-uniqueness of DMET fixed point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we are going to show that not only the projection but also the fixed point found in the DMET procedure is arbitrary and that it can be arbitrarily far away from the ''exact result''. In our case the ''exact result'' is the three-particle ground state of the following ''target'' Hamiltonian $\\boldsymbol{h}_{tar}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_tar=np.zeros((6,6))\n",
    "phi_tar_un=np.zeros((6))\n",
    "\n",
    "\n",
    "\n",
    "phi_tar[2:6,0]=copy.deepcopy(phi_env[:,0])\n",
    "\n",
    "phi_tar[:,1]=copy.deepcopy(phi[:,0])\n",
    "\n",
    "phi_tar[:,2]=copy.deepcopy(phi[:,1])\n",
    "\n",
    "phi_tar[:,3]=copy.deepcopy(phi[:,2])\n",
    "\n",
    "phi_tar[:,4]=copy.deepcopy(phi[:,3]-np.dot(phi_tar[:,0].transpose(),phi[:,3])*phi_tar[:,0])/ \\\n",
    "copy.deepcopy((np.dot(phi[:,3]-np.dot(phi_tar[:,0].transpose(),phi[:,3])*phi_tar[:,0], \\\n",
    "       phi[:,3]-np.dot(phi_tar[:,0].transpose(),phi[:,3])*phi_tar[:,0]))**0.5)\n",
    "phi_tar_un[:] = copy.deepcopy(phi[:,4])- np.dot(phi_tar[:,0].transpose(),phi[:,4])*phi_tar[:,0]-\\\n",
    "                np.dot(phi_tar[:,4].transpose(),phi[:,4])*phi_tar[:,4]\n",
    "phi_tar[:,5]=copy.deepcopy(phi_tar_un[:])/copy.deepcopy(np.dot(phi_tar_un,phi_tar_un))**0.5\n",
    "    \n",
    "\n",
    "epsilon_tar   = np.zeros((6))\n",
    "epsilon_tar[0]= -2\n",
    "epsilon_tar[1]= -1\n",
    "epsilon_tar[2]= -0.5\n",
    "epsilon_tar[3]=  0.5\n",
    "epsilon_tar[4]=  1\n",
    "epsilon_tar[5]=  2\n",
    "\n",
    "h_tar=np.zeros((6,6))\n",
    "h_tar=multi_dot([phi_tar,np.diag(epsilon_tar),phi_tar.transpose()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " where we have defined the orthogonal set of eigenfunctions as\n",
    " \n",
    " \\begin{align}\n",
    " \\{ \\phi_1^{\\rm tar}=\\varphi_1^{B},\n",
    "     \\phi_2^{\\rm tar}=\\phi_1,\n",
    "     \\phi_3^{\\rm tar}=\\phi_2,\n",
    "     \\phi_4^{\\rm tar}=\\phi_{3},\n",
    "     \\phi_5^{\\rm tar}=\\frac{1}{|\\phi_4-\\langle \\phi_1^{\\rm tar}|\\phi_4 \\rangle \\phi_1^{\\rm tar}|}(\\phi_4-\\langle \\phi_1^{\\rm tar}|\\phi_4\\rangle \\phi_1^{\\rm tar}),\\nonumber\\\\\n",
    "      \\phi_6^{\\rm tar}=\\frac{1}{|\\phi_5-\\langle \\phi_1^{\\rm tar}|\\phi_5 \\rangle \\phi_1^{\\rm tar}|}(\\phi_5-\\langle \\phi_1^{\\rm tar}|\\phi_5\\rangle \\phi_1^{\\rm tar})+\n",
    "      \\frac{1}{|\\phi_5-\\langle \\phi_5^{\\rm tar}|\\phi_5 \\rangle \\phi_5^{\\rm tar}|}(\\phi_5-\\langle \\phi_5^{\\rm tar}|\\phi_5\\rangle \\phi_5^{\\rm tar}\n",
    "     )  \\}    \n",
    " \\end{align}\n",
    "and $\\{\\phi_1,...,\\phi_5 \\}$ are the five lowest eigenstates of $\\boldsymbol{h_s}$. Further we have chosen\n",
    "\\begin{align*}\n",
    "    \\epsilon_1^{\\rm tar}=-2, \\epsilon_2^{\\rm tar}=-1, \\epsilon_3^{\\rm tar}=-0.5, \\epsilon_4^{\\rm tar}=0.5, \\epsilon_5^{\\rm tar}=1, \\epsilon_6^{\\rm tar}=2.\n",
    "\\end{align*}\n",
    "For this Hamiltonian the three particle ground state is $E^{tar}=-3.5$ and the 1RDM is \n",
    "$\\gamma^{tar}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.49796  0.35483  0.02354 -0.16463 -0.0215   0.3098 ]\n",
      " [ 0.35483  0.27354  0.08537 -0.0215   0.0598   0.24796]\n",
      " [ 0.02354  0.08537  0.3285   0.1098   0.44796  0.00483]\n",
      " [-0.16463 -0.0215   0.1098   0.89796 -0.04517  0.22354]\n",
      " [-0.0215   0.0598   0.44796 -0.04517  0.67354 -0.11463]\n",
      " [ 0.3098   0.24796  0.00483  0.22354 -0.11463  0.3285 ]]\n"
     ]
    }
   ],
   "source": [
    "gamma_tar=np.dot(phi_tar[:,0:3],phi_tar[:,0:3].transpose())\n",
    "print(gamma_tar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The impurity 1RDM submatrix reads:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.49796 0.35483]\n",
      " [0.35483 0.27354]]\n"
     ]
    }
   ],
   "source": [
    "gamma_imp=gamma_tar[0:2,0:2]\n",
    "print(gamma_imp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we assume that a DMET iteration step led us to an auxiliary Hamiltonian of the form of the Hamiltonian $\\boldsymbol{h}^{s}$ and determine the CAS orbitals like we did before we define the embedded Hamiltonian:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{h}^{'\\rm tar} =[\\mathbf{C}^{\\rm CAS}]^T \\mathbf{h}^{\\rm tar} \\mathbf{C}^{\\rm CAS}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.00000e+00 -1.00000e+00 -6.32456e-01  6.32456e-01]\n",
      " [-1.00000e+00  0.00000e+00 -6.32456e-01 -6.32456e-01]\n",
      " [-6.32456e-01 -6.32456e-01  6.00000e-01  4.44089e-16]\n",
      " [ 6.32456e-01 -6.32456e-01  4.44089e-16 -6.00000e-01]]\n"
     ]
    }
   ],
   "source": [
    "h_tar_prime=multi_dot([phi_CAS.transpose(),h,phi_CAS])\n",
    "print(h_tar_prime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Diagonalizing this Hamiltonian and keeping only the two lowest (embedded) orbitals we obtain an embedded 1RDM of the target Hamiltonian in the CAS basis as\n",
    "\\begin{align}\n",
    "    \\boldsymbol{\\gamma}'^{\\rm tar}_{\\rm CAS}&\\equiv\\sum_{k=1}^2 \\varphi_{k}^{\\rm emb}(\\mu)\\varphi_{k}^{\\rm emb}(\\nu)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.5      0.33333  0.26352 -0.26352]\n",
      " [ 0.33333  0.5      0.26352  0.26352]\n",
      " [ 0.26352  0.26352  0.16667 -0.     ]\n",
      " [-0.26352  0.26352 -0.       0.83333]]\n"
     ]
    }
   ],
   "source": [
    "n_emb, phi_emb=LA.eigh(h_tar_prime)\n",
    "gamma_tar_emb = np.dot(phi_emb[:,0:2],phi_emb[:,0:2].transpose()) \n",
    "print(np.around(gamma_tar_emb,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transforming the 1RDM into the site basis by\n",
    "\\begin{align}\n",
    "\\gamma'^{\\rm tar}_{\\rm emb}(i,j)=\\sum_{\\mu,\\nu=1}^4 \\gamma'^{\\rm tar}_{\\rm CAS}(\\mu,\\nu)\n",
    "\\varphi^{\\rm CAS}_{\\mu}(i)\\varphi^{\\rm CAS}_{\\nu}(j)\n",
    "\\end{align}\n",
    "leads to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.5      0.33333 -0.      -0.16667  0.       0.33333]\n",
      " [ 0.33333  0.5      0.33333 -0.      -0.16667 -0.     ]\n",
      " [-0.       0.33333  0.4      0.13333 -0.2     -0.26667]\n",
      " [-0.16667 -0.       0.13333  0.1     -0.06667 -0.2    ]\n",
      " [ 0.      -0.16667 -0.2     -0.06667  0.1      0.13333]\n",
      " [ 0.33333 -0.      -0.26667 -0.2      0.13333  0.4    ]]\n"
     ]
    }
   ],
   "source": [
    "gamma_tar_ij = multi_dot([phi_CAS,gamma_tar_emb,phi_CAS.transpose()])\n",
    "print(np.around(gamma_tar_ij,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We notice that the 1RDM $\\boldsymbol{\\gamma'}^{\\rm tar}_{emb}$ constructed from the embedded Hamiltonian does not agree with the target 1RDM $\\boldsymbol{\\gamma}^{\\rm tar}$ even on the impurity $A$. That is, the approximate impurity 1RDM is\n",
    "\\begin{align}\n",
    " \\boldsymbol{\\gamma}'^{\\rm tar}_{\\rm imp}=\\begin{pmatrix}\n",
    " 0.5 &0.33333 \\\\\n",
    " 0.33333 &0.5 \n",
    "\\end{pmatrix},\n",
    "\\end{align}\n",
    "while it does not agree with the exact impurity 1RDM we have reached self-consistency of the DMET algorithm since\n",
    "\\begin{align}\n",
    "    \\boldsymbol{\\gamma}^s_{\\rm imp} = \\boldsymbol{\\gamma}'^{\\rm tar}_{\\rm imp},\n",
    "\\end{align}\n",
    "and thus our DMET iteration is finished. <br> Additionally to finding a 1RDM that does not correspond to the actual 1RDM of the system, also the energy estimate is not necessarily good. To demonstrate this we are going to use the following formula to calculate first the energy of the fragment A:\n",
    "\\begin{align}\n",
    "\\label{eq:tot_energy_new}\n",
    "\\epsilon_{f}^{exact}=\\sum_{i=1,2,j=1-6}\\boldsymbol{h}^{tar}_{i,j}\\gamma^{tar}(j,i)=-0.55242\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.55242\n"
     ]
    }
   ],
   "source": [
    "epsilon_f=0\n",
    "for ii in range(2):\n",
    "    for jj in range(6):\n",
    "        epsilon_f=epsilon_f+h_tar[ii,jj]*gamma_tar[jj,ii]\n",
    "print(np.around(epsilon_f,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, because in practice we do not have the correct 1RDM that corresponds to this Hamiltonian available we need to calculate the fragment energy using the embedded 1RDM:\n",
    "\\begin{align}\n",
    "\\label{eq:tot_energy_appr}\n",
    "\\epsilon_{f}^{emb}=\\sum_{i=1,2,j=1-6}\\boldsymbol{h}^{tar}_{i,j}\\gamma^{'tar}_{emb}(j,i)=-0.43816\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.43816\n"
     ]
    }
   ],
   "source": [
    "epsilon_f_emb=0\n",
    "for ii in range(2):\n",
    "    for jj in range(6):\n",
    "        epsilon_f_emb=epsilon_f_emb+h_tar[ii,jj]*gamma_tar_ij[jj,ii]\n",
    "print(np.around(epsilon_f_emb,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the same procedure after adding to embedded 1RDM the environment orbital that we had originally discarded (so as to have three particles)\n",
    "\n",
    "\\begin{align}\n",
    "    \\gamma_{emb}^{tar,tot}(j,i)= \\boldsymbol{\\gamma}'^{\\rm tar}_{emb}(j,i)+\\phi_4^{*B}(j)\\phi^B_4(i)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_tot = np.zeros((6,6))\n",
    "gamma_tot = copy.deepcopy(gamma_tar_ij)\n",
    "gamma_tot[2:6,2:6] = copy.deepcopy(gamma_tot[2:6,2:6])+ np.dot(phi_env[:,3:4],phi_env[:,3:4].transpose())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we obtain the same wrong fragment energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.43816\n"
     ]
    }
   ],
   "source": [
    "epsilon_f_emb_tot=0\n",
    "for ii in range(2):\n",
    "    for jj in range(6):\n",
    "        epsilon_f_emb_tot=epsilon_f_emb_tot+h_tar[ii,jj]*gamma_tot[jj,ii]\n",
    "print(np.around(epsilon_f_emb_tot,5))"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}