{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "73d4113e",
   "metadata": {},
   "source": [
    "# Couple Cluster theory\n",
    "\n",
    "## Ground state\n",
    "\n",
    "The CC wavefunction can be written as:\n",
    "$$\n",
    "    \\Psi_{\\text{CC}} = e^T |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "The cluster operator $T$ can be expressed as:\n",
    "$$\n",
    "    T = T_1 + T_2 + T_3 + \\dots\n",
    "$$\n",
    "\n",
    "where the $T_n$ operator is defined as:\n",
    "\n",
    "$$\n",
    "    T_n = \\frac{1}{n!^2}\\sum t_{ijk\\dots}^{abc\\dots} a^\\dagger ib^\\dagger jc^\\dagger k\\dots\n",
    "$$\n",
    "\n",
    "For example, the single excitation operator $T_1$ is defined as:\n",
    "\n",
    "$$\n",
    "    T_1 = \\sum_{ia} t_i^a a^\\dagger i\n",
    "$$\n",
    "\n",
    "and the double excitation operator $T_2$ is defined as:\n",
    "\n",
    "$$\n",
    "    T_2 = \\frac{1}{4}\\sum_{ijab} t_{ij}^{ab} a^\\dagger ib^\\dagger j\n",
    "$$\n",
    "\n",
    "The exponential operator $e^T$ can be expanded as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    e^T  = & 1 + T + \\frac{1}{2}T^2 + \\frac{1}{3!}T^3 + \\dots \\\\\n",
    "      = & 1 + T_1 + T_2 + \\dots \\\\\n",
    "     & + \\frac{1}{2}T_1^2 + T_1T_2 + \\frac{1}{2}T_2^2 + \\dots \\\\\n",
    "     & + \\frac{1}{3!}T_1^3 + \\frac{1}{2}T_1^2T_2 + \\frac{1}{2}T_1T_2^2 + \\frac{1}{3!}T_2^3 + \\dots\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The single excitation can be expressed as:\n",
    "\n",
    "$$\n",
    "    \\Phi_1 = T_1 |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "and the double excitation can be expressed as:\n",
    "\n",
    "$$\n",
    "    \\Phi_2 = (T_2 + \\frac{1}{2}T_1^2 )|\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "Return to the Schrödinger equation, we can obtain that:\n",
    "\n",
    "$$\n",
    "    \\hat{H} \\Psi_\\text{CC} = \\hat{H} e^T |\\Phi_0\\rangle = E e^T |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "Multiplying the equation by $e^{-T}$ from the left, we can obtain:\n",
    "\n",
    "$$\n",
    "    e^{-T} \\hat{H} e^T |\\Phi_0\\rangle = E |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "We introduce the effective Hamiltonian (similarity transformed Hamiltonian):\n",
    "\n",
    "$$\n",
    "    \\bar{H} = e^{-T}H e^T\n",
    "$$\n",
    "\n",
    "Therefore, the equation can be rewritten as:\n",
    "\n",
    "$$\n",
    "    \\bar{H} |\\Phi_0\\rangle = E |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "By projecting onto excitation determinants, we can obtain:\n",
    "\n",
    "$$\n",
    "    \\langle \\Phi_{ij\\dots}^{ab\\dots} | \\bar{H} | \\Phi_0 \\rangle = 0\n",
    "$$\n",
    "\n",
    "where $\\Phi_{ij\\dots}^{ab\\dots}$ is the excitation determinant.\n",
    "\n",
    "The ground state energy can be expressed as:\n",
    "\n",
    "$$\n",
    "    E_{\\text{CC}} = \\langle \\Phi_0 | \\bar{H} | \\Phi_0 \\rangle\n",
    "$$\n",
    "\n",
    "The electronic Hamiltonian operator can be expressed as:\n",
    "$$\n",
    "    \\hat{H} = E_{\\text{HF}} +  \\sum_{pq} F_{pq} p^\\dagger q + \\frac{1}{4} \\sum_{pqrs} \\langle pq||rs \\rangle p^\\dagger q^\\dagger s r\n",
    "$$\n",
    "\n",
    "There are four terms in the one electron operator: OO, OV, VO, and VV. For the two-electron operator, there are 9 distinct terms: OOOO, OOOV, OVOO, OOVV, VVOO, OVOV, VOVV, VVVO, VVVV.\n",
    "\n",
    "The effective Hamiltonian can be rewritten using Baker-Campbell-Hausdorff formula as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\bar{H} = & \\hat{H} + \\left[\\hat{H}, T\\right] + \\frac{1}{2!}\\left[\\left[\\hat{H}, T\\right], T\\right] + \\dots  \\\\\n",
    "    = & \\hat{H} + \\left[\\hat{H}, T_1\\right] + \\left[\\hat{H}, T_2\\right] + \\frac{1}{2}\\left[\\left[\\hat{H}, T_1\\right], T_1\\right] \\\\\n",
    "    & + \\frac{1}{2}\\left[\\left[\\hat{H}, T_2\\right], T_2\\right] + \\left[\\left[\\hat{H}, T_1\\right], T_2\\right] + \\dots\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "knowing that:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\left[\\left[\\hat{H}, T_1\\right], T_2\\right]  = \\left[\\left[\\hat{H}, T_2\\right], T_1\\right]\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Only the single and double cluster amplitudes contribute to the CC energy (that means only consider the contribution of single and double), and therefore the energy can be expressed as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "     E_{\\text{CC}} & =  \\langle \\Phi_0 | \\bar{H} |\\Phi_0 \\rangle \\\\\n",
    "    & = \\langle \\Phi_0 | \\hat{H} | \\Phi_0 \\rangle + \\sum_{ia} F_{ia} t_i^a + \\frac{1}{2} \\sum_{ijab} \\langle ij||ab \\rangle t_{i}^{a} t_{j}^{b} + \\frac{1}{4} \\sum_{ijab} \\langle ij||ab \\rangle t_{ij}^{ab} \\\\\n",
    "    & = E_{\\text{HF}} + \\sum_{ia} F_{ia} t_i^a + \\frac{1}{4} \\sum_{ijab} \\langle ij||ab \\rangle(t_{ij}^{ab}+t_{i}^{a} t_{j}^{b} -t_{i}^{b} t_{j}^{a} )\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $F_{ia}$ is the Fock matrix element, and $\\langle ij||ab \\rangle$ is the two-electron integral.\n",
    "The $\\langle \\Phi_0 | \\hat{H} | \\Phi_0 \\rangle$ equals to the Hartree-Fock energy $E_{\\text{HF}}$.\n",
    "\n",
    "### CCSD\n",
    "\n",
    "In CCSD, the cluster operator is truncated at the single and double excitation level:\n",
    "\n",
    "$$\n",
    "    T = T_1 + T_2\n",
    "$$\n",
    "\n",
    "and the wavefunction can be expressed as:\n",
    "\n",
    "$$\n",
    "    |\\Psi_{\\text{CCSD}} \\rangle = e^{T_1+T_2} |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "Projecting onto singly and doubly excited determinants, we can obtain:\n",
    "\n",
    "$$\n",
    "    \\langle \\Phi_{i}^{a} | e^{-T_1}\\hat{H}e^{T_1} + \\left[e^{-T_1}\\hat{H}e^{T_1}, T_2\\right] | \\Phi_0 \\rangle = 0\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    & \\langle \\Phi_{ij}^{ab} | e^{-T_2-T_1}\\hat{H}e^{T_1+T_2} | \\Phi_0 \\rangle \n",
    "    = \\langle \\Phi_{ij}^{ab} | e^{-T_2}e^{-T_1}\\hat{H}e^{T_1}e^{T_2} | \\Phi_0 \\rangle \n",
    "    \\\\\n",
    "    & = \\langle \\Phi_{ij}^{ab} | e^{-T_1}\\hat{H}e^{T_1} + \\left[e^{-T_1}\\hat{H}e^{T_1}, T_2\\right] + \\frac{1}{2}\\left[\\left[e^{-T_1}\\hat{H}e^{T_1}, T_2\\right], T_2\\right] | \\Phi_0 \\rangle = 0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The energy can be expressed as:\n",
    "\n",
    "$$\n",
    "    E_{\\text{CCSD}} = \\sum_a\\sum_i t_i^a f_{ia} + \\frac{1}{4} \\sum_{ab}\\sum_{ij}\\tau_{ij}^{ab} \\langle ij||ab \\rangle\n",
    "$$\n",
    "\n",
    "### CC2\n",
    "\n",
    "In the CC2 method, the projection onto the singly excited determinant\n",
    "is the same as the CCSD method:\n",
    "$$\n",
    "    \\langle \\Phi_{i}^{a} | e^{-T_1}\\hat{H}e^{T_1} + \\left[e^{-T_1}\\hat{H}e^{T_1}, T_2\\right] | \\Phi_0 \\rangle = 0\n",
    "$$\n",
    "\n",
    "and the projection onto the doubly excited determinant is:\n",
    "\n",
    "$$\n",
    "    \\langle \\Phi_{ij}^{ab} | e^{-T_1}\\hat{H}e^{T_1} + \\left[F, T_2\\right] | \\Phi_0 \\rangle = 0\n",
    "$$\n",
    "\n",
    "## EOM-CC\n",
    "\n",
    "### EOM-EE (Electronically Excited States)\n",
    "\n",
    "Excited state can be expressed as:\n",
    "$$\n",
    "    |\\Psi_{\\text{ex}}\\rangle = \\mathbf{R} |\\Phi_\\text{g}\\rangle\n",
    "$$\n",
    "where $|\\Psi_{\\text{ex}}\\rangle$ is the excited state, and $|\\Phi_\\text{g}\\rangle$ is the ground state.\n",
    "$\\mathbf{R}$ is the excitation operator, which can be expressed as:\n",
    "$$\n",
    "    \\mathbf{R} = \\sum_{n} \\mathbf{R}_n\n",
    "$$\n",
    "\n",
    "$$\n",
    "    \\mathbf{R}_n = \\frac{1}{n!^2}\\sum r_{ijk\\dots}^{abc\\dots} a^\\dagger ib^\\dagger jc^\\dagger k\\dots\n",
    "$$\n",
    "\n",
    "The ground state is approximated by:\n",
    "\n",
    "$$\n",
    "    |\\Phi_\\text{g}\\rangle = e^T |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "where $|\\Phi_0\\rangle$ is an arbitrary Slater determinant, usually SCF solution.\n",
    "\n",
    "Therefore, the excited state can be expressed as:\n",
    "\n",
    "$$\n",
    "    |\\Psi_{\\text{ex}}\\rangle =\\mathbf{R} e^T  |\\Phi_0\\rangle = e^T \\mathbf{R} |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "which means the right excitation operator $\\mathbf{R}$ and the cluster operator $T$ are commutative.\n",
    "\n",
    "\n",
    "Inserting the equation into the Schrödinger equation, we can obtain:\n",
    "\n",
    "$$\n",
    "    H \\mathbf{R} e^T  |\\Phi_0\\rangle = E_\\text{ex} \\mathbf{R} e^T  |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "The $T$ operator and $\\mathbf{R}$ operator are necessarily commutative, so we can rewrite the equation\n",
    "using the effective Hamiltonian operator $\\bar{H}= e^{-T}H e^T$:\n",
    "\n",
    "$$\n",
    "    e^{-T}H e^T \\mathbf{R}  |\\Phi_0\\rangle = \\bar{H} \\mathbf{R} |\\Phi_0\\rangle =\n",
    "     E_\\text{ex}  \\mathbf{R}  |\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "and then we can obtain the equation:\n",
    "\n",
    "$$\n",
    "   (\\bar{H}-E_\\text{ex}) \\mathbf{R}  |\\Phi_0\\rangle = 0\n",
    "$$\n",
    "\n",
    "\n",
    "The goal of any EOM-CC calculation is to determine the energy difference between the initial and target states:\n",
    "\n",
    "$$\n",
    "\\omega_k = E_k - E_{\\text{CC}}\n",
    "$$\n",
    "\n",
    "With the commutation relation between the cluster operator $T$ and the excitation operator $\\mathbf{R}$, we can obtain the EOM-CC equation:\n",
    "\n",
    "$$\n",
    "[\\bar{H},\\mathbf{R}]|\\Phi_0\\rangle = \\omega_k \\mathbf{R}|\\Phi_0\\rangle\n",
    "$$\n",
    "\n",
    "The effective Hamiltonian can be expressed in spin-orbital basis as:\n",
    "\n",
    "$$\n",
    "    \\bar{H} = \\sum_{pq} \\mathbf{F}_{pq} a^\\dagger_p a_q + \\frac{1}{4} \\sum_{pqrs} \\mathbf{W}_{pqsr} a^\\dagger_p a^\\dagger_q a_s a_r\n",
    "$$\n",
    "\n",
    "Note that the CC wave operator $e^T$ is not unitary, so the effective Hamiltonian $\\bar{H}$ is not Hermitian. As a result, each root of $\\bar{H}$ is associated with two eigenvectors, which are the left and right eigenvectors. The left eigenvector is defined as:\n",
    "\n",
    "$$\n",
    "    \\langle \\tilde{\\Psi} | = \\langle \\Psi_L | =\n",
    "    \\langle \\Phi_0 | \\mathbf{L} e^{-T}\n",
    "$$\n",
    "\n",
    "and the right eigenvector is defined as:\n",
    "\n",
    "$$\n",
    "    | \\Psi \\rangle = | \\Psi_R \\rangle =\n",
    "    e^T \\mathbf{R} |\\Phi_0 \\rangle\n",
    "$$\n",
    "\n",
    "Note that the left deexcitation operator $\\mathbf{L}$ is not commutative with the cluster operator $T$. The left eigenfunction is defined as:\n",
    "\n",
    "$$\n",
    "\\langle \\Psi_0 | \\mathbf{L}\\bar{H} = \\langle \\Phi_0 | \\mathbf{L}\\omega_k\n",
    "$$\n",
    "\n",
    "The left de-excitation operator $\\mathbf{L}$ can be expressed as:\n",
    "\n",
    "$$\n",
    "    \\mathbf{L} = \\sum_{n} \\mathbf{L}_n\n",
    "$$\n",
    "\n",
    "$$\n",
    "    \\mathbf{L}_n = \\frac{1}{n!^2}\\sum l_{abc\\dots}^{ijk\\dots} i^\\dagger aj^\\dagger bk^\\dagger c\\dots\n",
    "$$\n",
    "\n",
    "The left and right operators satisfy the biorthogonality condition:\n",
    "\n",
    "$$\n",
    "    \\langle \\mathbf{L}_i | \\mathbf{R}_j \\rangle = \\delta_{ij}\n",
    "$$\n",
    "\n",
    "The two sets of solutions satisfy the biorthogonality condition:\n",
    "\n",
    "$$\n",
    "    \\langle \\tilde{\\Psi}^{(i)} | \\Psi^{(j)} \\rangle = 1\n",
    "$$\n",
    "\n",
    "If the C is unity, we can rewrite the equation as:\n",
    "\n",
    "$$\n",
    "    \\langle \\tilde{\\Psi} | \\Psi \\rangle = 1\n",
    "$$\n",
    "\n",
    "and therefore the energy can be expressed as:\n",
    "\n",
    "$$\n",
    "    E = \\langle \\tilde{\\Psi} | H | \\Psi \\rangle =\n",
    "    \\langle \\Phi_0 |\\mathbf{L} \\bar{H} \\mathbf{R}| \\Phi_0 \\rangle\n",
    "$$\n",
    "\n",
    "Note that for ground state, we have $\\mathbf{L} = \\mathbf{R} = 1$, so the energy can be expressed as:\n",
    "\n",
    "$$\n",
    "    E_{\\text{CC}} = \\langle \\Phi_0 | \\bar{H} | \\Phi_0 \\rangle\n",
    "$$\n",
    "\n",
    "The eigenvalues can be obtained by diagonalizing the matrix $\\mathbf{A}$ in the basis of the reference, single excited and doubly excited determinants. The Jacobian matrix can be expressed as:\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \n",
    "\\begin{pmatrix}\n",
    "\\langle\\Phi_i^a|[\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger k\\} ]| \\Phi_0\\rangle &  \\langle\\Phi_i^a|[\\tilde{H},\\{c^\\dagger d^\\dagger kl\\} ]| \\Phi_0\\rangle \\\\\n",
    "\\langle\\Phi_{ij}^{ab}|[\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger k\\}] | \\Phi_0\\rangle &  \\langle\\Phi_{ij}^{ab}|\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger d^\\dagger kl\\} | \\Phi_0\\rangle\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "For CC2, the Jacobian matrix can be expressed as:\n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \n",
    "\\begin{pmatrix}\n",
    "\\langle\\Phi_i^a|[\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger k\\} ]| \\Phi_0\\rangle &  \\langle\\Phi_i^a|[\\tilde{H},\\{c^\\dagger d^\\dagger kl\\} ]| \\Phi_0\\rangle \\\\\n",
    "\\langle\\Phi_{ij}^{ab}|[\\tilde{H},\\{c^\\dagger k\\}] | \\Phi_0\\rangle &  \\langle\\Phi_{ij}^{ab}|F,\\{c^\\dagger d^\\dagger kl\\} | \\Phi_0\\rangle\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "And the whole eigenfunction can be expressed as:\n",
    "\n",
    "$$\n",
    "\\mathbf{A} \n",
    "\\begin{pmatrix}\n",
    "R_1 \\\\ R_2\n",
    "\\end{pmatrix}\n",
    "= \\omega \n",
    "\\begin{pmatrix}\n",
    "R_1 \\\\ R_2\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "L_1 & L_2\n",
    "\\end{pmatrix}\n",
    "\\mathbf{A} \n",
    "= \\begin{pmatrix}\n",
    "L_1 & L_2\n",
    "\\end{pmatrix}\\omega \n",
    "$$\n",
    "\n",
    "Then, the right and left eigenfunctions can be expressed as:\n",
    "\n",
    "$$\n",
    "\\sum_{kc}\\langle\\Phi_i^a|[\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger k\\} ]| \\Phi_0\\rangle r_k^c + \\sum_{klcd} \\langle\\Phi_i^a|[\\tilde{H},\\{c^\\dagger d^\\dagger kl\\} ]| \\Phi_0\\rangle r_{kl}^{cd} = \\omega r_i^a\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\sum_{kc}\\langle\\Phi_{ij}^{ab}|[\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger k\\}] | \\Phi_0\\rangle r_k^c + \\sum_{klcd} \\langle\\Phi_{ij}^{ab}|\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger d^\\dagger kl\\} | \\Phi_0\\rangle = \\omega r_{ij}^{ab}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\sum_{ia}l_a^i \\langle\\Phi_i^a|[\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger k\\} ]| \\Phi_0\\rangle + \\sum_{ijab} l_{ij}^{ab} \n",
    "\\langle\\Phi_{ij}^{ab}|[\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger k\\}] | \\Phi_0\\rangle = \\omega l_k^c\n",
    "$$\n",
    "\n",
    "$$\n",
    " \\sum_{ia}l_a^i\\langle\\Phi_i^a|[\\tilde{H},\\{c^\\dagger d^\\dagger kl\\} ]| \\Phi_0\\rangle + \\sum_{ijab} l_{ij}^{ab} \n",
    "\\langle\\Phi_{ij}^{ab}|\\tilde{H}+[\\tilde{H},T_2],\\{c^\\dagger d^\\dagger kl\\} | \\Phi_0\\rangle = \\omega l_{kl}^{cd}\n",
    "$$\n",
    "\n",
    "\n",
    "### EOM-EA (Electron Attachment)\n",
    "\n",
    "For EOM-EA state, the right excitation operator $\\mathbf{R}$ is defined as:\n",
    "\n",
    "$$\n",
    "    \\mathbf{R}^{\\text{EA}} = \\sum_{a} r_a a^\\dagger + \\frac{1}{2} \\sum_{iab} r_{i}^{ab} ia^\\dagger b^\\dagger + \\cdots\n",
    "$$\n",
    "\n",
    "\n",
    "### EOM-IP (Ionization Potential)\n",
    "\n",
    "For EOM-IP state, the right excitation operator $\\mathbf{R}$ is defined as:\n",
    "\n",
    "$$\n",
    "    \\mathbf{R}^{\\text{IP}} = \\sum_{i} r_i i + \\frac{1}{2} \\sum_{ija} r_{ij}^{a} ija^\\dagger + \\cdots\n",
    "$$\n",
    "\n",
    "\n",
    "### EOM-SF (Spin-Flip)\n",
    "\n",
    "For EOM-SF state, the right excitation operator $\\mathbf{R}$ is defined as:\n",
    "\n",
    "$$\n",
    "    \\mathbf{R}_{M_s=-1} \\equiv \\mathbf{R}^{\\text{IP}} = \\sum_{ia}r_i^aa_\\beta^\\dagger i_\\alpha + \\frac{1}{2}\\sum_{ijab}r_{ij}^{ab}a_\\beta^\\dagger i_\\alpha b_\\sigma^\\dagger j_\\sigma + \\cdots\n",
    "$$\n",
    "\n",
    "where $\\sigma=\\alpha$ or $\\beta$.\n",
    "\n",
    "## Derivative\n",
    "\n",
    "### General\n",
    "\n",
    "The asymmetrized, perturbation-independent, deexcitation operator\n",
    "$\\Lambda$ can be expressed as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& \\Lambda = \\sum_{n} \\Lambda_n \\\\\n",
    "& \\Lambda_1 = \\sum_{ia} \\lambda_i^a i^\\dagger a \\\\\n",
    "& \\Lambda_2 = \\frac{1}{4}\\sum_{ijab} \\lambda_{ij}^{ab} a_i^\\dagger a_j^\\dagger a_b a_a \\\\\n",
    "& \\Lambda_n = \\frac{1}{n!^2}\\sum \\lambda_{abc\\dots}^{ijk\\dots} i^\\dagger aj^\\dagger bk^\\dagger c\\dots\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The $\\lambda$ amplitude satisfies the condition:\n",
    "\n",
    "$$\n",
    "    \\langle \\Phi_0 | \\Lambda | \\Phi \\rangle\n",
    "    = - \\langle \\Phi_0 | (H_Ne^T)_c | \\Phi \\rangle\n",
    "    \\langle \\Phi | (H_Ne^T)_c - \\Delta E| \\Phi \\rangle ^{-1}\n",
    "$$\n",
    "\n",
    "Therefore, we have:\n",
    "\n",
    "$$\n",
    "    \\langle \\Phi_0 | (H_N^\\xi e^T)_c | \\Phi_0 \\rangle + \\langle \\Phi_0 | \\Lambda | \\Phi \\rangle \\langle \\Phi | (H_N^\\xi e^T)_c | \\Phi_0 \\rangle  = \\Delta E^\\xi\n",
    "$$\n",
    "\n",
    "Finally the derivative of the energy can be expressed as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\Delta E^\\xi = & \\mathbf{D}_{ij}f_{ij}^\\xi + \\mathbf{D}_{ab}f_{ab}^\\xi + \\sum_{pqrs}\\Gamma(pq,rs) \\langle pq||rs \\rangle^\\xi \\\\\n",
    "    = & \\mathbf{D}_{ij}f_{ij}^\\xi + \\mathbf{D}_{ab}f_{ab}^\\xi  \\\\\n",
    "   & + \\Gamma(ij,ab)\\langle ij||ab\\rangle^\\xi\n",
    "    + \\Gamma(ab,cd)\\langle ab||cd\\rangle^\\xi   \\\\\n",
    "   & + \\Gamma(kl,ij)\\langle kl||ij\\rangle^\\xi\n",
    "    + \\Gamma(ja,bi)\\langle ja||bi\\rangle^\\xi   \\\\\n",
    "   & + \\Gamma(ai,bc)\\langle ai||bc\\rangle^\\xi\n",
    "    + \\Gamma(jk,ia)\\langle jk||ia\\rangle^\\xi   \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where the $D_{pq}$ is the relaxed density matrix, and the $\\Gamma_{pq,rs}$ is the effective density matrix.\n",
    "\n",
    "### CCSD Derivative\n",
    "\n",
    "The CCSD gradient can be expressed as:\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial E_{\\text{CCSD}}}{\\partial \\xi} =\n",
    "\\sum_{pq} D_{pq}f_{pq}^{(\\xi)} + \\sum_{pqrs}\\Gamma(pq,rs) \\langle pq||rs \\rangle^\\xi\n",
    "+ \\sum_{pq}I_{pq}S_{pq}^{\\xi}\n",
    "$$\n",
    "\n",
    "### EOM-CC Derivative\n",
    "\n",
    "The derivative of the EOM-CC energy can be expressed as:\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial E_{\\text{EOM-CC}}}{\\partial \\xi} =\n",
    "D_{pq}f_{pq}^{(\\xi)} + \\frac{1}{4}\\Gamma_{pqrs} \\langle pq||rs \\rangle^\\xi\n",
    "+ I_{pq}S_{pq}^{\\xi}\n",
    "$$\n",
    "\n",
    "where:\n",
    "overlap derivative:\n",
    "$$\n",
    "S_{pq}^{\\xi} = c_{\\mu p} \\frac{\\partial S_{\\mu\\nu}}{\\partial \\xi} c_{\\nu q}\n",
    "= c_{\\mu p} c_{\\nu q} \\left( \\langle \\frac{\\partial\\chi_\\mu}{\\partial\\xi}|\\chi_\\nu \\rangle + \\langle \\chi_\\mu|\\frac{\\partial\\chi_\\nu}{\\partial\\xi} \\rangle \\right)\n",
    "$$\n",
    "\n",
    "ERI derivative:\n",
    "$$\n",
    "\\langle pq||rs \\rangle^\\xi = c_{\\mu p} c_{\\nu q} c_{\\lambda r} c_{\\sigma s} \\frac{\\partial \\langle \\mu\\lambda||\\nu\\sigma \\rangle}{\\partial \\xi}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\phi_p = c_{\\mu p} \\xi_{\\mu}\n",
    "$$\n",
    "\n",
    "Fock matrix derivative:\n",
    "$$\n",
    "f_{pq}^{(\\xi)} = c_{\\mu p}c_{\\nu q} \\left(\\frac{\\partial h_{\\mu\\nu}}{\\partial\\xi}+ c_{\\lambda i}c_{\\sigma i}\\frac{\\partial\\langle \\mu\\lambda||\\nu\\sigma \\rangle}{\\partial\\xi} \\right)\n",
    "= h_{pq}^\\xi + \\sum_{m} \\langle pm||qm \\rangle^\\xi\n",
    "$$\n",
    "\n",
    "Derivative of the MO coefficients:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial c_{\\mu p}}{\\partial \\xi} = \\sum_q U_{qp}^\\xi c_{\\mu q}\n",
    "$$\n",
    "\n",
    "where $U_{qp}^\\xi$ is the CPHF coefficient.\n",
    "\n",
    "Relaxed density matrix:\n",
    "$$\n",
    "D_{pq} = \\gamma_{pq} + z_{pq}\n",
    "$$\n",
    "\n",
    "Effective density matrix:\n",
    "\n",
    "$$\n",
    "\\gamma_{pq} = \\langle \\Phi_0 | \\mathbf{L}[p^\\dagger qe^T]_c\\mathbf{R} |\\Phi_0 \\rangle\n",
    "+ \\langle \\Phi_0 | \\mathcal{Z}[p^\\dagger qe^T]_c|\\Phi_0 \\rangle\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\Gamma_{pqrs}\n",
    "= \\langle \\Phi_0 | \\mathbf{L}[p^\\dagger q^\\dagger sr e^T]_c\\mathbf{R} |\\Phi_0 \\rangle\n",
    "+ \\langle \\Phi_0 | \\mathcal{Z}[p^\\dagger q^\\dagger sr e^T]_c|\\Phi_0 \\rangle\n",
    "$$\n",
    "\n",
    "where $c$ indicates the limitation to connected diagrams.\n",
    "\n",
    "Auxiliary deexcitation operator $\\mathcal{Z}$:\n",
    "\n",
    "$$\n",
    "\\mathcal{Z} = \\sum_{n} \\mathcal{Z}_n\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathcal{Z}_n = \\frac{1}{n!^2}\\sum \\zeta_{abc\\dots}^{ijk\\dots} i^\\dagger aj^\\dagger bk^\\dagger c\\dots\n",
    "$$\n",
    "\n",
    "### Lagrangian of EOM-CC derivative\n",
    "\n",
    "The EOM energy can be expressed as:\n",
    "\n",
    "\n",
    "\n",
    "The full energy derivative can be expressed as:\n",
    "\n",
    "$$\n",
    "\\frac{\\mathrm{d} E }{\\mathrm{d} \\xi} = \\frac{\\partial E}{\\partial \\xi} +\n",
    "\\frac{\\partial E}{\\partial L}\\frac{\\partial L}{\\partial \\xi} + \\frac{\\partial E}{\\partial R}\\frac{\\partial R}{\\partial \\xi} +\n",
    "\\frac{\\partial E}{\\partial T}\\frac{\\partial T}{\\partial \\xi} + \\frac{\\partial E}{\\partial C}\\frac{\\partial C}{\\partial \\xi}\n",
    "$$\n",
    "\n",
    "The first term is the Hellmann-Feynman contribution:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& \\frac{\\partial E}{\\partial \\xi} =  \\sum_{pq} h_{pq}^\\xi \\gamma_{pq}' + \\frac{1}{4}\\sum_{pqrs}\\langle pq||rs \\rangle^\\xi \\Gamma_{pqrs}' \n",
    "\\\\\n",
    "& h_{pq}^\\xi = \\frac{\\partial h_{pq}}{\\partial \\xi} = \\sum_{\\mu\\nu}C_{\\mu p}h_{\\mu\\nu}^\\xi C_{\\nu q} \n",
    "\\\\\n",
    "& h_{\\mu\\nu}^\\xi = \\langle \\chi_\\mu | \\frac{\\partial \\hat{h}}{\\partial \\xi} | \\chi_\\nu \\rangle + \\langle \\frac{\\partial \\chi_\\mu}{\\partial\\xi} | \\hat{h} | \\chi_\\nu \\rangle + \\langle \\chi_\\mu | \\hat{h} | \\frac{\\partial \\chi_\\nu}{\\partial\\xi} \\rangle\n",
    "\\\\\n",
    "& \\langle pq||rs \\rangle^\\xi = \\frac{\\partial \\langle pq||rs \\rangle}{\\partial \\xi} = \\sum_{\\mu\\nu\\lambda\\sigma}C_{\\mu p}C_{\\nu q}\\langle \\chi_\\mu\\chi_\\nu||\\chi_\\lambda\\chi_\\sigma \\rangle^\\xi C_{\\lambda r}C_{\\sigma s} \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\langle \\chi_\\mu\\chi_\\nu||\\chi_\\lambda\\chi_\\sigma \\rangle^\\xi = &\n",
    "\\langle \\frac{\\partial \\chi_\\mu}{\\partial\\xi}\\chi_\\nu||\\chi_\\lambda\\chi_\\sigma \\rangle +\n",
    "\\langle \\chi_\\mu\\frac{\\partial \\chi_\\nu}{\\partial\\xi}||\\chi_\\lambda\\chi_\\sigma \\rangle  \n",
    "\\\\\n",
    "& + \\langle \\chi_\\mu\\chi_\\nu||\\frac{\\partial \\chi_\\lambda}{\\partial\\xi}\\chi_\\sigma \\rangle +\n",
    "\\langle \\chi_\\mu\\chi_\\nu||\\chi_\\lambda\\frac{\\partial \\chi_\\sigma}{\\partial\\xi} \\rangle\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The EOM energy is stationary with respect to the left and right eigenvectors, \n",
    "so the second and third terms are zero:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial E}{\\partial L} = \\frac{\\partial E}{\\partial R} = 0\n",
    "$$\n",
    "\n",
    "Then there is so-called amplitude response $\\frac{\\partial E}{\\partial T}$ \n",
    " and orbital response $\\frac{\\partial E}{\\partial C}$.\n",
    "\n",
    "The Lagrangian derivative can be expressed as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\partial \\mathcal{L}(L, R, T, C,Z, \\Lambda, \\Omega)}{\\partial \\xi} = & \n",
    "\\langle \\Phi_0 Le^{-T} | \\frac{\\partial H}{\\partial\\xi} | e^TR\\Phi_0 \\rangle \n",
    "\\\\\n",
    "+ & \\langle \\Phi_0 Ze^{-T} | \\frac{\\partial H}{\\partial\\xi} | e^T\\Phi_0 \\rangle \n",
    "\\\\ \n",
    "+ & \\frac{1}{2}\\sum_{pq}\\lambda_{pq}\\frac{\\partial f_{pq}}{\\partial\\xi} + \n",
    "\\sum_{pq}\\omega_{pq}\\frac{\\partial S_{pq}}{\\partial\\xi}\n",
    "\\\\\n",
    "= & \\sum_{pq}h_{pq}^\\xi\\rho_{pq} + \\frac{1}{4}\\sum_{pqrs}\\langle pq||rs \\rangle^\\xi \\Pi_{pqrs} \\\\\n",
    "+ & \\sum_{pq}\\omega_{pq}S_{pq}^{\\xi}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The effective density matrices $\\rho$ and $\\Pi$ can be expressed as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& \\rho = \\gamma' + \\gamma'' + \\gamma'''\n",
    "\\\\\n",
    "& \\Pi = \\Gamma' + \\Gamma'' + \\Gamma'''\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where the $\\gamma'$ and $\\Gamma'$ are the so-called non-relaxed density matrices:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& \\gamma' = \\frac{1}{2}\\langle \\Psi_L | p^\\dagger q+q^\\dagger p |\\Psi_R \\rangle\n",
    "\\\\\n",
    "& \\Gamma' = \\frac{1}{2}\\langle \\Psi_L | p^\\dagger q^\\dagger sr + s^\\dagger r^\\dagger pq |\\Psi_R \\rangle\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "and $\\gamma''$ and $\\Gamma''$ are amplitude response contributions:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& \\gamma'' = \\frac{1}{2}\\langle \\Phi_0Ze^{-T} | p^\\dagger q+q^\\dagger p |e^T\\Phi_0 \\rangle\n",
    "\\\\\n",
    "& \\Gamma'' = \\frac{1}{2}\\langle \\Phi_0Ze^{-T} | p^\\dagger q^\\dagger sr + s^\\dagger r^\\dagger pq |e^T\\Phi_0 \\rangle\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "and $\\gamma'''$ and $\\Gamma'''$ are orbital response contributions. $\\gamma'''$\n",
    "is related to the Lagrange multiplier $\\lambda$, and $\\Gamma'''$ is related to \n",
    "$\\gamma'''$ and $\\delta$.\n",
    "\n",
    "## Properties\n",
    "\n",
    "Different CC2 variants:\n",
    "\n",
    "![IP](figures/IP.png)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}