{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test of F12 matrix elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cd \"/home/scemama/QCaml\";;\n",
    "#use \"topfind\";;\n",
    "#require \"jupyter.notebook\";;\n",
    "\n",
    "let png_image name = \n",
    "    Jupyter_notebook.display_file ~base64:true \"image/png\" (\"Notebooks/images/\"^name)\n",
    ";;\n",
    "\n",
    "#require \"lacaml.top\";;\n",
    "#require \"alcotest\";;\n",
    "#require \"str\";;\n",
    "#require \"bigarray\";;\n",
    "#require \"zarith\";;\n",
    "#require \"getopt\";;\n",
    "#directory \"_build\";;\n",
    "#directory \"_build/Basis\";;\n",
    "#directory \"_build/CI\";;\n",
    "#directory \"_build/MOBasis\";;\n",
    "#directory \"_build/Nuclei\";;\n",
    "#directory \"_build/Parallel\";;\n",
    "#directory \"_build/Perturbation\";;\n",
    "#directory \"_build/SCF\";;\n",
    "#directory \"_build/Utils\";;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modules to load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load \"Constants.cmo\";;\n",
    "#load_rec \"Util.cma\";;\n",
    "#load_rec \"Matrix.cmo\";;\n",
    "#load_rec \"Simulation.cmo\";;\n",
    "#load_rec \"Spindeterminant.cmo\";;\n",
    "#load_rec \"Determinant.cmo\";;\n",
    "#load_rec \"HartreeFock.cmo\";;\n",
    "#load_rec \"MOBasis.cmo\";;\n",
    "#load_rec \"F12CI.cmo\";;\n",
    "\n",
    "open Lacaml.D;;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Printers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#install_printer AngularMomentum.pp_string ;;\n",
    "#install_printer Basis.pp ;;\n",
    "#install_printer Charge.pp ;;\n",
    "#install_printer Coordinate.pp ;;\n",
    "#install_printer Vector.pp;;\n",
    "#install_printer Matrix.pp;;\n",
    "#install_printer Util.pp_float_2darray;;\n",
    "#install_printer Util.pp_float_array;;\n",
    "#install_printer Util.pp_matrix;;\n",
    "#install_printer HartreeFock.pp ;;\n",
    "#install_printer Fock.pp ;;\n",
    "#install_printer MOClass.pp ;;\n",
    "let pp_mo ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;\n",
    "#install_printer pp_mo;;\n",
    "#install_printer DeterminantSpace.pp;;\n",
    "#install_printer SpindeterminantSpace.pp;;\n",
    "#install_printer Bitstring.pp;;\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let basis_filename  = \"/home/scemama/qp2/data/basis/6-31g\" \n",
    "let aux_basis_filename = \"/home/scemama/qp2/data/basis/cc-pvdz\" \n",
    "let nuclei = Nuclei.of_zmt_string \"be\" \n",
    "let frozen_core = false\n",
    "let multiplicity = 1\n",
    "let state = 1\n",
    "\n",
    "let basis =  Basis.of_nuclei_and_basis_filenames ~nuclei [basis_filename] \n",
    "let aux_basis =  Basis.of_nuclei_and_basis_filenames ~nuclei (basis_filename :: aux_basis_filename :: []) \n",
    "let f12 = F12factor.gaussian_geminal 1.0 \n",
    "let charge = 0 \n",
    "\n",
    "\n",
    "let simulation =\n",
    "    Simulation.make \n",
    "        ~f12 ~charge ~multiplicity ~nuclei\n",
    "        ~cartesian:true\n",
    "        basis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hartree-Fock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let hf =  HartreeFock.make ~guess:`Hcore  simulation  ;;\n",
    "\n",
    "let mo_basis = MOBasis.of_hartree_fock hf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FCI-F12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notations:\n",
    "\n",
    "* $\\langle ij || kl \\rangle = \\int \\phi_i(r_1) \\phi_j(r_2) \\frac{1}{r_{12}} \\phi_k(r1) \\phi_l(r2) $ \n",
    "* $\\left[ ij || kl \\right] = \\int \\phi_i(r_1) \\phi_j(r_2) f_{12} \\phi_k(r1) \\phi_l(r2) $ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Common functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let f12 = Util.of_some @@ Simulation.f12 simulation \n",
    "\n",
    "let mo_num = MOBasis.size mo_basis \n",
    "\n",
    "let pp_spindet = Spindeterminant.pp mo_num\n",
    "\n",
    "let pp_det = Determinant.pp mo_num\n",
    "\n",
    ";;\n",
    "\n",
    "#install_printer pp_spindet ;;\n",
    "#install_printer pp_det ;;\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let simulation_aux = \n",
    "    let charge        = Charge.to_int @@ Simulation.charge simulation \n",
    "    and multiplicity  = Electrons.multiplicity @@ Simulation.electrons simulation\n",
    "    and nuclei        = Simulation.nuclei simulation\n",
    "    in\n",
    "    let general_basis = \n",
    "      Basis.general_basis @@ Simulation.basis simulation\n",
    "    in\n",
    "    GeneralBasis.combine [\n",
    "      general_basis ; GeneralBasis.read aux_basis_filename\n",
    "    ]\n",
    "    |> Basis.of_nuclei_and_general_basis nuclei\n",
    "    |> Simulation.make ~f12 ~charge ~multiplicity ~nuclei \n",
    "\n",
    "\n",
    "let aux_basis = \n",
    "  MOBasis.of_mo_basis simulation_aux mo_basis\n",
    "\n",
    "let aux_num = \n",
    "  MOBasis.size aux_basis\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let () = ignore @@ MOBasis.f12_ints   aux_basis\n",
    "let () = ignore @@ MOBasis.two_e_ints aux_basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integral-based functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\langle I | \\hat{H} | J \\rangle = \\begin{cases}\n",
    "\\sum_i h_{ii} + \\frac{1}{2} \\sum_{ij} \\langle ij || ij \\rangle \\text{ when } |J\\rangle = |I\\rangle \\\\\n",
    "h_{ik} + \\sum_{j} \\langle ij || kj \\rangle \\text{ when } |J\\rangle = \\hat{T}_i^k |I\\rangle \\\\\n",
    "\\langle ij || kl \\rangle \\text{ when } |J\\rangle = \\hat{T}_{ij}^{kl} |I\\rangle \\\\\n",
    "\\end{cases}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle I | \\hat{F} | J \\rangle = \\begin{cases}\n",
    "\\sum_i f_{ii} + \\frac{1}{2} \\sum_{ij} \\langle ij || ij \\rangle \\text{ when } |J\\rangle = |I\\rangle \\\\\n",
    "f_{ik} + \\sum_{j} \\langle ij || kj \\rangle \\text{ when } |J\\rangle = \\hat{T}_i^k |I\\rangle \\\\\n",
    "\\langle ij || kl \\rangle \\text{ when } |J\\rangle = \\hat{T}_{ij}^{kl} |I\\rangle \\\\\n",
    "\\end{cases}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let cancel_singles = false \n",
    "\n",
    "let mos_cabs = \n",
    "  Util.list_range (mo_num+1) aux_num\n",
    "  \n",
    "let mos_abs = \n",
    "  Util.list_range 1 aux_num\n",
    "  \n",
    "let mos_in = \n",
    "  Util.list_range 1 mo_num\n",
    "\n",
    "let mos_a k =\n",
    "  Determinant.alfa k\n",
    "  |> Spindeterminant.to_list\n",
    "  \n",
    "let mos_b k =\n",
    "  Determinant.beta k\n",
    "  |> Spindeterminant.to_list\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### H integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let h_one =\n",
    "  let h = \n",
    "    MOBasis.one_e_ints aux_basis\n",
    "  in fun i j _ -> h.{i,j}\n",
    "  \n",
    "let h_two = \n",
    "  let two_e_ints = MOBasis.two_e_ints aux_basis in\n",
    "  let h2 i j k l (s:Spin.t) (s':Spin.t) =\n",
    "    if s' <> s then\n",
    "        ERI.get_phys two_e_ints i j k l\n",
    "    else\n",
    "        (ERI.get_phys two_e_ints i j k l) -.                                                 \n",
    "        (ERI.get_phys two_e_ints i j l k)\n",
    "  in\n",
    "  h2\n",
    "\n",
    "    \n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### F12 integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let f_two = \n",
    "  let two_e_ints = MOBasis.f12_ints aux_basis in\n",
    "  let f2 i j k l (s:Spin.t) (s':Spin.t) =\n",
    "    if s' <> s then\n",
    "        0.5 *. F12.get_phys two_e_ints i j k l\n",
    "    else\n",
    "        0.5 *. (\n",
    "         (F12.get_phys two_e_ints i j k l) -.                                                 \n",
    "         (F12.get_phys two_e_ints i j l k) )\n",
    "  in\n",
    "  let f3 i j k l (s:Spin.t) (s':Spin.t) = \n",
    "      if (i=k && j<>l) || (j=l && i<>k) then\n",
    "        0.\n",
    "      else\n",
    "        f2 i j k l s s'\n",
    "  in\n",
    "  if cancel_singles then f3 else f2\n",
    "\n",
    "let f_one = fun _ _ _ -> 0.\n",
    "\n",
    "(*\n",
    "let f_two = h_two\n",
    "\n",
    "let f_one = h_one\n",
    "*)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determinant-based functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let f12_integrals mo_basis =\n",
    "  ( f_one, f_two, None )\n",
    "\n",
    "let h_ij mo_basis ki kj =\n",
    "  let integrals =\n",
    "    List.map (fun f -> f mo_basis)\n",
    "      [ CI.h_integrals ]\n",
    "  in\n",
    "  CIMatrixElement.make integrals ki kj \n",
    "  |> List.hd\n",
    "\n",
    "\n",
    "let f_ij mo_basis ki kj =\n",
    "  let integrals =\n",
    "    List.map (fun f -> f mo_basis)\n",
    "      [ f12_integrals ]\n",
    "  in\n",
    "  CIMatrixElement.make integrals ki kj \n",
    "  |> List.hd                                                                                  \n",
    "\n",
    "\n",
    "let hf_ij mo_basis ki kj =\n",
    "  let integrals =\n",
    "    List.map (fun f -> f mo_basis)\n",
    "      [ CI.h_integrals ; f12_integrals ]\n",
    "  in\n",
    "  CIMatrixElement.make integrals ki kj\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determinant space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let is_a_double det_space =\n",
    "    let mo_class = DeterminantSpace.mo_class det_space in\n",
    "    let mo_num = Array.length @@ MOClass.mo_class_array mo_class in\n",
    "    let m l =\n",
    "      List.fold_left (fun accu i ->\n",
    "          let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)\n",
    "        ) (Bitstring.zero mo_num) l\n",
    "    in\n",
    "    let aux_mask   = m (MOClass.auxiliary_mos mo_class) in\n",
    "    fun k ->                                                                                  \n",
    "      let alfa =\n",
    "        Determinant.alfa k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let beta =\n",
    "        Determinant.beta k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let a = Bitstring.logand aux_mask alfa\n",
    "      and b = Bitstring.logand aux_mask beta\n",
    "      in\n",
    "      match Bitstring.popcount a + Bitstring.popcount b with\n",
    "      | 2 | 1 -> true\n",
    "      | 0 | _ -> false\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let in_space = \n",
    "    DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num\n",
    "    \n",
    "let aux_space = \n",
    "    DeterminantSpace.fci_of_mo_basis aux_basis ~frozen_core \n",
    "\n",
    "let det_space_in () =\n",
    "    DeterminantSpace.determinant_stream in_space\n",
    "\n",
    "let det_space_out () =\n",
    "    let s = \n",
    "        DeterminantSpace.determinant_stream aux_space\n",
    "    in\n",
    "    Stream.from (fun _ ->\n",
    "        try\n",
    "          let is_a_double = is_a_double in_space in\n",
    "          let rec result () =\n",
    "            let ki = Stream.next s in\n",
    "            if is_a_double ki then\n",
    "              Some (ki,ki)\n",
    "            else\n",
    "              result ()\n",
    "          in\n",
    "          result ()\n",
    "        with Stream.Failure -> None\n",
    "      )\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let ci = CI.make ~n_states:state in_space\n",
    "\n",
    "let ci_coef, ci_energy = Lazy.force ci.eigensystem \n",
    "\n",
    "let _ = print_newline () \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Permutation operator $p_{12}$ that generates a new determinant with electrons 1 and 2 swapped."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "let p12 det_space =\n",
    "    let mo_class = DeterminantSpace.mo_class det_space in\n",
    "    let mo_num = Array.length @@ MOClass.mo_class_array mo_class in\n",
    "    let m l =\n",
    "      List.fold_left (fun accu i ->\n",
    "          let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)\n",
    "        ) (Bitstring.zero mo_num) l\n",
    "    in\n",
    "    let aux_mask = m (MOClass.auxiliary_mos mo_class) in\n",
    "    let not_aux_mask =\n",
    "      Bitstring.(shift_left_one mo_num (mo_num-1) |> minus_one |> logxor aux_mask)\n",
    "    in\n",
    "    fun k ->\n",
    "      let alfa =\n",
    "        Determinant.alfa k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let beta =\n",
    "        Determinant.beta k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let a = Bitstring.logand aux_mask alfa\n",
    "      and b = Bitstring.logand aux_mask beta\n",
    "      in\n",
    "      match Bitstring.popcount a, Bitstring.popcount b with\n",
    "      | 2, 0  \n",
    "      | 0, 2 -> Some (Determinant.negate_phase k)\n",
    "      | 1, 1 -> Some (Determinant.of_spindeterminants\n",
    "                  (Spindeterminant.of_bitstring  @@\n",
    "                    Bitstring.(logor b (logand not_aux_mask alfa)) )\n",
    "                  (Spindeterminant.of_bitstring  @@\n",
    "                    Bitstring.(logor a (logand not_aux_mask beta))\n",
    "                  ) )\n",
    "      | 1, 0 \n",
    "      | 0, 1 -> Some k\n",
    "      | _ -> None\n",
    "                    \n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrices $\\langle I | H | \\alpha \\rangle$ and $\\langle I | F | \\alpha \\rangle$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let out_list =\n",
    "    Util.stream_to_list (det_space_out ())\n",
    "    \n",
    "let in_list =\n",
    "    Util.stream_to_list (det_space_in ())\n",
    "    \n",
    "let det_a = Array.of_list out_list\n",
    "            |> Array.map (fun (i,_) -> i)\n",
    "\n",
    "let det_I = Array.of_list in_list\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "let m_H_aux, m_F_aux = \n",
    "        \n",
    "    let rec col_vecs_list accu_H accu_F = function\n",
    "    | [] -> \n",
    "        List.rev accu_H, \n",
    "        List.rev accu_F \n",
    "    | (ki, ki') :: rest when ki = ki' ->\n",
    "        begin\n",
    "            let h, f = \n",
    "              List.map (fun kj ->\n",
    "              match hf_ij aux_basis kj ki with\n",
    "              | [ a ; b ] -> a, b\n",
    "              | _ -> assert false ) in_list\n",
    "              |> List.split\n",
    "            in\n",
    "            let h  = Vec.of_list h in\n",
    "            let f  = Vec.of_list f in\n",
    "            scal 0.5 f;  (* Specific to He singlet *)\n",
    "            col_vecs_list (h::accu_H) (f::accu_F) rest\n",
    "        end\n",
    "    | (ki, ki') :: rest ->\n",
    "        begin\n",
    "        assert false;\n",
    "        end\n",
    "    in\n",
    "    let h, f = \n",
    "      col_vecs_list [] [] out_list\n",
    "    in\n",
    "    Mat.of_col_vecs_list h,\n",
    "    Mat.of_col_vecs_list f\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to generate all intermediate external determinants $|\\alpha \\rangle$ between $|I\\rangle$ and $|J\\rangle$, with a positive phase factor.\n",
    "\n",
    "* `exc` is the degree of excitation between $|I\\rangle$ and $|J\\rangle$\n",
    "* `cabs` is the number of electrons in the CABS\n",
    "* `l` is the degree of excitation between $|I\\rangle$ and $|\\alpha \\rangle$\n",
    "* `r` is the degree of excitation between $|\\alpha \\rangle$ and $|J\\rangle$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let generate_alphas ki kj exc cabs l r =\n",
    "  (* Check input excitation degree *)\n",
    "  let d = Determinant.degree ki kj in\n",
    "  if d <> exc then\n",
    "    Printf.sprintf \"Invalid excitation degree. Expected %d and found %d.\" exc d\n",
    "    |> failwith;\n",
    "    \n",
    "  (* Generate single excitations *)\n",
    "  let all_singles ki =\n",
    "      let mos_a, mos_b = Determinant.to_lists ki in\n",
    "      [ List.map (fun hole -> \n",
    "          List.map (fun particle ->\n",
    "                  Determinant.single_excitation Spin.Alfa hole particle ki\n",
    "              ) mos_abs\n",
    "          ) mos_a\n",
    "      ;\n",
    "        List.map (fun hole -> \n",
    "          List.map (fun particle ->\n",
    "                  Determinant.single_excitation Spin.Beta hole particle ki\n",
    "              ) mos_abs\n",
    "          ) mos_b \n",
    "      ]\n",
    "      |> List.concat\n",
    "      |> List.concat\n",
    "      |> List.filter (fun x -> Determinant.is_none x = false)\n",
    "      |> List.map (fun x -> Determinant.set_phase Phase.Pos x)\n",
    "  in\n",
    "  \n",
    "  (* Generate double excitations *)\n",
    "  let all_doubles ki =\n",
    "      let mos_a, mos_b = Determinant.to_lists ki in\n",
    "      [ List.map (fun hole ->  (* Alpha-Alpha *)\n",
    "          List.map (fun particle ->\n",
    "            List.map (fun hole' -> \n",
    "              List.map (fun particle' ->\n",
    "                if hole' > hole && particle' < particle then\n",
    "                  Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Alfa hole' particle' ki)\n",
    "                else None\n",
    "                      ) mos_abs\n",
    "                  ) mos_a \n",
    "              ) mos_abs\n",
    "          ) mos_a\n",
    "        ; \n",
    "        List.map (fun hole ->  (* Beta-Beta *)\n",
    "          List.map (fun particle ->\n",
    "            List.map (fun hole' -> \n",
    "              List.map (fun particle' ->\n",
    "                if hole' > hole && particle' < particle then\n",
    "                  Some (Determinant.double_excitation Spin.Beta hole particle Spin.Beta hole' particle' ki)\n",
    "                else None\n",
    "                      ) mos_abs\n",
    "                  ) mos_b \n",
    "              ) mos_abs\n",
    "          ) mos_b\n",
    "        ;\n",
    "        List.map (fun hole ->  (* Alpha-Beta *)\n",
    "          List.map (fun particle ->\n",
    "            List.map (fun hole' -> \n",
    "              List.map (fun particle' ->\n",
    "                  Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Beta hole' particle' ki)\n",
    "                      ) mos_abs\n",
    "                  ) mos_b \n",
    "              ) mos_abs\n",
    "          ) mos_a\n",
    "      ]\n",
    "      |> List.concat\n",
    "      |> List.concat\n",
    "      |> List.concat\n",
    "      |> List.concat\n",
    "      |> Util.list_some\n",
    "      |> List.filter (fun x -> Determinant.is_none x = false)\n",
    "      |> List.map (fun x -> Determinant.set_phase Phase.Pos x)\n",
    "  in\n",
    "\n",
    "  (* Generate left and right excitations *)\n",
    "  let al = \n",
    "      match l with\n",
    "      | 1 -> all_singles ki\n",
    "      | 2 -> all_doubles ki\n",
    "      | _ -> assert false\n",
    "  in\n",
    "  let ar = \n",
    "      match r with\n",
    "      | 1 -> all_singles kj\n",
    "      | 2 -> all_doubles kj\n",
    "      | _ -> assert false\n",
    "  in\n",
    "  \n",
    "  let mo_class = DeterminantSpace.mo_class in_space in\n",
    "  let m l =\n",
    "    List.fold_left (fun accu i ->\n",
    "        let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)\n",
    "      ) (Bitstring.zero mo_num) l\n",
    "  in\n",
    "  let aux_mask = m (MOClass.auxiliary_mos mo_class) in\n",
    "  let good_cabs k = \n",
    "      let alfa =\n",
    "        Determinant.alfa k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let beta =\n",
    "        Determinant.beta k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let a = Bitstring.logand aux_mask alfa\n",
    "      and b = Bitstring.logand aux_mask beta\n",
    "      in\n",
    "      Bitstring.popcount a + Bitstring.popcount b = cabs\n",
    "  in\n",
    "  let good_lr k =\n",
    "    Determinant.degree ki k = l &&\n",
    "    Determinant.degree k kj = r \n",
    "  in\n",
    "\n",
    "\n",
    "  \n",
    "  (* Merge lists in a set of unique determinants *)\n",
    "  List.concat [ al; ar ]\n",
    "  |> List.sort_uniq compare\n",
    "  \n",
    "  (* Filter out all determinants with incorrect numbers of electrons in the CABS *)\n",
    "  |> List.filter good_cabs \n",
    "  \n",
    "  (* Filter out all determinants with incorrect excitation numbers *)\n",
    "  |> List.filter good_lr \n",
    "  \n",
    "  \n",
    "let compute_HaaF ki alphas kj =\n",
    "   List.fold_left (fun accu alpha -> accu\n",
    "      +. h_ij aux_basis ki alpha\n",
    "      *. f_ij aux_basis alpha kj\n",
    "    ) 0. alphas\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let check n integral_value exc cabs lexc rexc =\n",
    "    let det_list =\n",
    "        match n with\n",
    "          | 0 -> det_I\n",
    "          | n -> Array.sub det_I 0 n\n",
    "    in\n",
    "    let result =\n",
    "        Printf.printf \"Checking ... \\n%!\";\n",
    "        let percent = ref 0 in\n",
    "        Array.mapi (fun i ki -> \n",
    "            let p = (10 * (i+1))/(Array.length det_list) in\n",
    "            if p <> !percent then\n",
    "                ( percent := p ; Printf.printf \" - %3d %%\\n%!\" (p*10) );\n",
    "            Array.mapi (fun j kj -> \n",
    "                if i > j then (0,0,0.,0.) else\n",
    "                if Determinant.degree ki kj <> exc then (0,0,0.,0.) else\n",
    "                let alphas = generate_alphas ki kj exc cabs lexc rexc in\n",
    "                let det_value =\n",
    "                  compute_HaaF ki alphas kj\n",
    "                in\n",
    "                let int_value = \n",
    "                  integral_value ki kj\n",
    "                in\n",
    "                (i,j,det_value,int_value)\n",
    "            ) det_list\n",
    "        ) det_list\n",
    "        |> Array.to_list\n",
    "        |> Array.concat\n",
    "  in\n",
    "  let i,j,d,v = \n",
    "      let rec aux k imax jmax emax dmax vmax = \n",
    "        if k = -1 then\n",
    "          imax, jmax, dmax, vmax\n",
    "        else\n",
    "          let i, j, d, v = result.(k) in\n",
    "          let e = abs_float (d -. v) in\n",
    "          if e > emax then\n",
    "            aux (k-1) i j e d v\n",
    "          else\n",
    "            aux (k-1) imax jmax emax dmax vmax\n",
    "      in\n",
    "      aux (n-1) 0 0 0. 0. 0.\n",
    "  in\n",
    "  let error = abs_float (d -. v) in\n",
    "  if error < epsilon_float then\n",
    "      Printf.printf \"OK: %e\\n%!\" error\n",
    "  else\n",
    "      Printf.printf \"Failed: (%d, %d) | %e %e | %e\\n%!\" i j d v error\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "******"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. 0 1 11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"0_1_11.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$i,j$ : orbital indices of MOs occupied in $|I\\rangle$.\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_{a} \\left( \n",
    "        \\sum_{i} \\sum_{j}\n",
    "            \\left(h_{ia} +  \\langle i j || a j \\rangle \\right) \n",
    "            \\left(f_{ai} +   \\left[ a j || i j \\right] \\right) + \n",
    "        \\sum_{i} \\sum_{\\bar{j}}\n",
    "            \\left(h_{ia} +  \\langle i \\bar{j} | a \\bar{j} \\rangle \\right) \n",
    "            \\left(f_{ai} +   \\left[ a \\bar{j} | i \\bar{j} \\right] \\right) \n",
    "    \\right) + \\\\\n",
    "\\sum_{\\bar{a}} \\left(\n",
    "        \\sum_{\\bar{i}} \\sum_{\\bar{j}}\n",
    "            \\left(h_{\\bar{i}\\bar{a}} +  \\langle \\bar{i} \\bar{j} || \\bar{a} \\bar{j} \\rangle \\right) \n",
    "            \\left(f_{\\bar{a}\\bar{i}} +   \\left[ \\bar{a} \\bar{j} || \\bar{i} \\bar{j} \\right] \\right) + \n",
    "        \\sum_{\\bar{i}} \\sum_{j}\n",
    "            \\left(h_{\\bar{i}\\bar{a}} +  \\langle \\bar{i} j | \\bar{a} j \\rangle \\right)\n",
    "            \\left(f_{\\bar{a}\\bar{i}} +   \\left[ \\bar{a} j | \\bar{i} j \\right] \\right) \n",
    "        \\right) \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "   (* Alpha-Beta *)\n",
    "   let s, s' = Spin.(Alfa, Beta) in \n",
    "   List.fold_left (fun accu j -> accu +. \n",
    "       List.fold_left (fun accu i -> accu +. \n",
    "           List.fold_left (fun accu a ->\n",
    "               accu +.( (h_one i a s  +. h_two i j a j s s') *. \n",
    "                        (f_one a i s  +. f_two a j i j s s') )\n",
    "                    +.( (h_one j a s' +. h_two i j i a s s') *. \n",
    "                        (f_one a j s' +. f_two i a i j s s') )\n",
    "             ) 0. mos_cabs\n",
    "         ) 0. (mos_a ki)\n",
    "     ) 0. (mos_b ki)\n",
    "   +.\n",
    "   (* Alpha-Alpha / Beta-Beta *)\n",
    "   List.fold_left (fun accu (mos,s) -> \n",
    "     let s' = s in accu +. \n",
    "     List.fold_left (fun accu j -> accu +. \n",
    "         List.fold_left (fun accu i -> accu +. \n",
    "               List.fold_left (fun accu a -> \n",
    "                   accu +.( (h_one i a s +. h_two i j a j s s') *. \n",
    "                            (f_one a i s +. f_two a j i j s s') )\n",
    "                   ) 0. mos_cabs\n",
    "           ) 0. mos\n",
    "       ) 0. mos\n",
    "     ) 0. [ (mos_a ki,Spin.Alfa) ; (mos_b ki,Spin.Beta)]\n",
    "     \n",
    "let _ =\n",
    "  check 0 integral_value 0 1 1 1\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. 0 1 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"0_1_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$i,j$ : orbital indices of MOs   occupied in $|I\\rangle$.\n",
    "\n",
    "$k$   : orbital indices of MOs unoccupied in $|I\\rangle$.\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_{a} \\sum_{k} \\sum_{i} \\sum_{j<i}\n",
    "            \\langle i j || a k \\rangle \n",
    "             \\left[ a k || i j \\right] + \n",
    "\\sum_{a} \\sum_{\\bar{k}} \\sum_{i} \\sum_{\\bar{j}}\n",
    "            \\langle i \\bar{j} | a \\bar{k} \\rangle \n",
    "             \\left[ a \\bar{k} | i \\bar{j} \\right]  + \\\\\n",
    "\\sum_{\\bar{a}} \\sum_{k} \\sum_{\\bar{i}} \\sum_{j}\n",
    "            \\langle \\bar{i} j | \\bar{a} k \\rangle \n",
    "             \\left[ \\bar{a} k | \\bar{i} j \\right] + \n",
    "\\sum_{\\bar{a}} \\sum_{\\bar{k}} \\sum_{\\bar{i}} \\sum_{\\bar{j}<\\bar{i}}\n",
    "            \\langle \\bar{i} \\bar{j} || \\bar{a} \\bar{k} \\rangle \n",
    "             \\left[ \\bar{a} \\bar{k} || \\bar{i} \\bar{j} \\right] \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "   (* mos unoccupied in both I and J *)\n",
    "   let mos_virt_a, mos_virt_b = \n",
    "     Array.init mo_num (fun i -> Some (i+1)) , \n",
    "     Array.init mo_num (fun i -> Some (i+1)) \n",
    "   in\n",
    "   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);\n",
    "   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);\n",
    "   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);\n",
    "   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);\n",
    "   \n",
    "   let mos_virt_a, mos_virt_b = \n",
    "     Array.to_list mos_virt_a |> Util.list_some,\n",
    "     Array.to_list mos_virt_b |> Util.list_some\n",
    "   in\n",
    "   \n",
    "   (* Alpha-Beta *)\n",
    "   let s, s' = Spin.(Alfa, Beta) in\n",
    "   List.fold_left (fun accu j -> accu +. \n",
    "       List.fold_left (fun accu i -> accu +. \n",
    "           List.fold_left (fun accu k -> accu +. \n",
    "               List.fold_left (fun accu a ->\n",
    "                   accu +. h_two i j a k s s'  *. f_two a k i j s s'\n",
    "                 ) 0. mos_cabs\n",
    "             ) 0. mos_virt_b\n",
    "         ) 0. (mos_a ki)\n",
    "     ) 0. (mos_b ki)\n",
    "   +.\n",
    "   List.fold_left (fun accu j -> accu +. \n",
    "       List.fold_left (fun accu i -> accu +. \n",
    "           List.fold_left (fun accu k -> accu +. \n",
    "               List.fold_left (fun accu a ->\n",
    "                   accu +. h_two i j a k s s'  *. f_two a k i j s s'\n",
    "                 ) 0. mos_cabs\n",
    "             ) 0. mos_virt_a\n",
    "         ) 0. (mos_b ki)\n",
    "     ) 0. (mos_a ki)\n",
    "   +.\n",
    "   (* Alpha-Alpha / Beta-Beta *)\n",
    "   List.fold_left (fun accu (mos_virt,mos,s) -> \n",
    "     let s' = s in accu +.\n",
    "     List.fold_left (fun accu j -> accu +. \n",
    "         List.fold_left (fun accu i -> accu +. \n",
    "             if i < j then accu else\n",
    "                 List.fold_left (fun accu k -> accu +. \n",
    "                     List.fold_left (fun accu a -> \n",
    "                           accu +. h_two i j a k s s'  *. f_two a k i j s s'\n",
    "                       ) 0. mos_cabs\n",
    "                   ) 0. mos_virt\n",
    "           ) 0. mos\n",
    "       ) 0. mos\n",
    "     ) 0. [ (mos_virt_a,mos_a ki,Spin.Alfa) ; (mos_virt_b,mos_b ki,Spin.Beta)]\n",
    "     \n",
    "let _ =\n",
    "  check 0 integral_value 0 1 2 2\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. 0 2 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "png_image \"0_2_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$i,j$ : orbital indices of MOs   occupied in $|I\\rangle$.\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_{a} \\sum_{b>a} \\sum_{i} \\sum_{j<i}\n",
    "            \\langle i j || a b \\rangle \n",
    "             \\left[ a b || i j \\right] + \n",
    "\\sum_{a} \\sum_{\\bar{b}} \\sum_{i} \\sum_{\\bar{j}}\n",
    "            \\langle i \\bar{j} | a \\bar{b} \\rangle  \n",
    "             \\left[ a \\bar{b} | i \\bar{j} \\right]  + \\\\\n",
    "\\sum_{\\bar{a}} \\sum_{b} \\sum_{\\bar{i}} \\sum_{j}\n",
    "            \\langle \\bar{i} j | \\bar{a} b \\rangle  \n",
    "             \\left[ \\bar{a} b | \\bar{i} j \\right] + \n",
    "\\sum_{\\bar{a}} \\sum_{\\bar{b}>\\bar{a}} \\sum_{\\bar{i}} \\sum_{\\bar{j}<\\bar{i}}\n",
    "            \\langle \\bar{i} \\bar{j} || \\bar{a} \\bar{b} \\rangle \n",
    "             \\left[ \\bar{a} \\bar{b} || \\bar{i} \\bar{j} \\right] \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "   (* Alpha-Beta *)\n",
    "   let s, s' = Spin.(Alfa, Beta) in\n",
    "   List.fold_left (fun accu j -> accu +. \n",
    "       List.fold_left (fun accu i -> accu +. \n",
    "           List.fold_left (fun accu b -> accu +. \n",
    "               List.fold_left (fun accu a ->\n",
    "                   accu +. h_two i j a b s s'  *. f_two a b i j s s'\n",
    "                 ) 0. mos_cabs\n",
    "             ) 0. mos_cabs\n",
    "         ) 0. (mos_a ki)\n",
    "     ) 0. (mos_b ki)\n",
    "   +.\n",
    "   (* Alpha-Alpha / Beta-Beta *)\n",
    "   List.fold_left (fun accu (mos,s) -> \n",
    "     let s' = s in accu +.\n",
    "     List.fold_left (fun accu j -> accu +. \n",
    "         List.fold_left (fun accu i -> accu +. \n",
    "             if i < j then accu else\n",
    "                 List.fold_left (fun accu b -> accu +. \n",
    "                     List.fold_left (fun accu a -> \n",
    "                         if b > a then accu else\n",
    "                           accu +. h_two i j a b s s'  *. f_two a b i j s s'\n",
    "                       ) 0. mos_cabs\n",
    "                   ) 0. mos_cabs\n",
    "           ) 0. mos\n",
    "       ) 0. mos\n",
    "     ) 0. [ (mos_a ki,Spin.Alfa) ; (mos_b ki,Spin.Beta)]\n",
    "     \n",
    "let _ =\n",
    "  check 100 integral_value 0 2 2 2\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. *1 1 11*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"1_1_11.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{i}^{k}|I\\rangle$\n",
    "\n",
    "$j$ : orbital indices of MOs occupied both in $|I\\rangle$  and $|J\\rangle$.\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_{a} \\sum_{j}\n",
    "            \\left( h_{ia} + \\langle i j || a j \\rangle \\right)\n",
    "            \\left( f_{ak} +  \\left[ a j || k j \\right] \\right) + \n",
    "\\sum_{\\bar{j}}\n",
    "            \\left( h_{ia} + \\langle i \\bar{j} | a \\bar{j} \\rangle \\right) \n",
    "            \\left( f_{ak} +  \\left[ a \\bar{j} | k \\bar{j} \\right] \\right)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "    let i, k, s, phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Single (phase, { hole ; particle ; spin })) ->\n",
    "            hole, particle, spin, phase\n",
    "          | _ -> assert false\n",
    "    in\n",
    "\n",
    "    let mos, mos' =\n",
    "      match s with\n",
    "      | Spin.Alfa -> mos_a ki, mos_b ki\n",
    "      | Spin.Beta -> mos_b ki, mos_a ki\n",
    "    in\n",
    "\n",
    "    let result = \n",
    "       (* Alpha-Beta *)\n",
    "       List.fold_left (fun accu j ->\n",
    "           let s' = Spin.other s in\n",
    "           accu +. List.fold_left (fun accu a ->\n",
    "                   accu +. (h_one i a s +. h_two i j a j s s')\n",
    "                        *. (f_one a k s +. f_two a j k j s s')\n",
    "                 ) 0. mos_cabs\n",
    "         ) 0. mos'\n",
    "       +.\n",
    "       (* Alpha-Alpha / Beta-Beta *)\n",
    "       List.fold_left (fun accu j ->\n",
    "           if j = i then accu else\n",
    "           let s' = s in\n",
    "           accu +. List.fold_left (fun accu a ->\n",
    "                   accu +. (h_one i a s +. h_two i j a j s s')\n",
    "                        *. (f_one a k s +. f_two a j k j s s')\n",
    "                 ) 0. mos_cabs\n",
    "         ) 0. mos\n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 50 integral_value 1 1 1 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let ki = det_I.(0)\n",
    "let kj = det_I.(0)\n",
    "\n",
    "\n",
    "let _ =\n",
    "    let alphas =\n",
    "        generate_alphas ki kj 0 1 2 2\n",
    "    in\n",
    "    compute_HaaF ki alphas kj\n",
    "\n",
    "let _ = integral_value ki kj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. *1 1 12*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"1_1_12.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. *1 1 21*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"1_1_21.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. *1 1 22*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"1_1_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. 1 2 22 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"1_2_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{i}^{k}|I\\rangle$\n",
    "\n",
    "$j$ : orbital indices of MOs occupied both in $|I\\rangle$  and $|J\\rangle$.\n",
    "\n",
    "\n",
    "$$\n",
    "\\sum_{j} \\sum_{b}\\sum_{a<b}\n",
    "    \\langle i j || a b \\rangle \n",
    "    \\left[ ab || kj \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "    let h, p, s, phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Single (phase, { hole ; particle ; spin })) ->\n",
    "            hole, particle, spin, phase\n",
    "          | _ -> assert false\n",
    "    in\n",
    "\n",
    "    let mos, mos' =\n",
    "      match s with\n",
    "      | Spin.Alfa -> mos_a ki, mos_b ki\n",
    "      | Spin.Beta -> mos_b ki, mos_a ki\n",
    "    in\n",
    "\n",
    "    let result = \n",
    "       (* Alpha-Beta *)\n",
    "       List.fold_left (fun accu j ->\n",
    "           let s' = Spin.other s in\n",
    "           accu +. List.fold_left (fun accu b ->\n",
    "               accu +. List.fold_left (fun accu a ->\n",
    "                   accu +. h_two h j a b s s'  *. f_two a b p j s s'\n",
    "                 ) 0. mos_cabs\n",
    "             ) 0. mos_cabs\n",
    "         ) 0. mos' \n",
    "       +.\n",
    "       (* Alpha-Alpha / Beta-Beta *)\n",
    "       List.fold_left (fun accu j ->\n",
    "           let s' = s in \n",
    "           accu +. List.fold_left (fun accu b ->\n",
    "             accu +. List.fold_left (fun accu a -> \n",
    "                 if b > a then accu else\n",
    "                   accu +. h_two h j a b s s'  *. f_two a b p j s s'\n",
    "               ) 0. mos_cabs\n",
    "           ) 0. mos_cabs\n",
    "         ) 0. mos \n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 50 integral_value 1 2 2 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9. *2 1 12*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"2_1_12.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10. *2 1 21*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"2_1_21.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 11. *2 1 22*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"2_1_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12. 2 2 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"2_2_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
    "\n",
    "$$\n",
    "\\sum_{b}\\sum_{a<b}\n",
    "    \\langle i j || a b \\rangle \n",
    "    \\left[ ab || kj \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "    let h, h', p, p', s, s', phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Double (phase,\n",
    "              { hole=h ; particle=p ; spin=s },\n",
    "              { hole=h'; particle=p'; spin=s'}) )  -> h, h', p, p', s, s', phase\n",
    "          | _ -> assert false\n",
    "    in\n",
    "\n",
    "    let result = \n",
    "       if s <> s' then (* Alpha-Beta *)\n",
    "           List.fold_left (fun accu b ->\n",
    "               accu +. List.fold_left (fun accu a ->\n",
    "                   accu +. h_two h h' a b s s'  *. f_two a b p p' s s'\n",
    "                 ) 0. mos_cabs\n",
    "             ) 0. mos_cabs\n",
    "       else (* Alpha-Alpha / Beta-Beta *)\n",
    "           List.fold_left (fun accu b ->\n",
    "             accu +. List.fold_left (fun accu a -> \n",
    "                 if b > a then accu else\n",
    "                   accu +. h_two h h' a b s s'  *. f_two a b p p' s s'\n",
    "               ) 0. mos_cabs\n",
    "           ) 0. mos_cabs\n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 100 integral_value 2 2 2 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 13. *3 1 2 2*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "png_image \"3_1_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{ijm}^{kln}|I\\rangle$\n",
    "\n",
    "$$\n",
    "\\sum_{a} \\langle i j || k a \\rangle  \\left[ m a || n l \\right] \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "    let h1, h2, h3, p1, p2, p3, s1, s2, s3, phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Triple (phase,\n",
    "              { hole=h1 ; particle=p1 ; spin=s1 },\n",
    "              { hole=h2 ; particle=p2 ; spin=s2 },\n",
    "              { hole=h3 ; particle=p3 ; spin=s3 }) )  -> h1, h2, h3, p1, p2, p3, s1, s2, s3, phase\n",
    "          | _ -> assert false\n",
    "    in\n",
    "\n",
    "(*\n",
    "123\n",
    "213\n",
    "231\n",
    "321\n",
    "312\n",
    "132\n",
    "*)\n",
    "\n",
    "    let result = \n",
    "       List.fold_left (fun accu a -> accu +.\n",
    "           h_two h1 h2 p1 a s1 s2  *. f_two h3 a p3 p2 s3 s2 -.\n",
    "           h_two h2 h1 p2 a s2 s1  *. f_two h3 a p3 p1 s3 s1 +.\n",
    "           h_two h2 h3 p2 a s2 s3  *. f_two h1 a p1 p3 s1 s3 -.\n",
    "           h_two h3 h2 p3 a s3 s2  *. f_two h1 a p1 p2 s1 s2 +.\n",
    "           h_two h3 h1 p3 a s3 s1  *. f_two h2 a p2 p1 s2 s1 -.\n",
    "           h_two h1 h3 p1 a s1 s3  *. f_two h2 a p2 p3 s2 s3 \n",
    "           ) 0. mos_cabs\n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 100 integral_value 3 1 2 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "OCaml default",
   "language": "OCaml",
   "name": "ocaml-jupyter"
  },
  "language_info": {
   "codemirror_mode": "text/x-ocaml",
   "file_extension": ".ml",
   "mimetype": "text/x-ocaml",
   "name": "OCaml",
   "nbconverter_exporter": null,
   "pygments_lexer": "OCaml",
   "version": "4.07.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}