QCaml/Notebooks/F12_matrix.ipynb

{
 "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": [
    "let png_image = print_endline ;;\n",
    "\n",
    "(* --------- *)\n",
    "\n",
    "#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"
   ]
  },
  {
   "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\";;"
   ]
  },
  {
   "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",
    "(*\n",
    "#install_printer DeterminantSpace.pp;;\n",
    "*)\n",
    "#install_printer SpindeterminantSpace.pp;;\n",
    "#install_printer Bitstring.pp;;\n",
    "\n",
    "(* --------- *)\n",
    "\n",
    "open Lacaml.D\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 \"c\" \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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let n_elec_alfa, n_elec_beta, n_elec = \n",
    "  let e =  Simulation.electrons simulation in\n",
    "  Electrons.(n_alfa e, n_beta e, n_elec e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hartree-Fock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let hf =  HartreeFock.make ~guess:`Hcore ~max_scf:10 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",
    "\n",
    "let f_one = h_one"
   ]
  },
  {
   "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": {
    "scrolled": false
   },
   "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": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "let randomize = false\n",
    "\n",
    "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 =\n",
    "  if randomize then\n",
    "    let n = 123456789 in\n",
    "    in_list\n",
    "    |> List.map (fun k -> (Random.int n, k))\n",
    "    |> List.sort compare\n",
    "    |> List.map  (fun (_,k) -> k)\n",
    "    |> Array.of_list\n",
    "  else\n",
    "    Array.of_list in_list\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/markdown"
   },
   "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"
   ]
  },
  {
   "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",
    "               if hole = particle then None else\n",
    "                  Some (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",
    "               if hole = particle then None else\n",
    "                 Some (Determinant.single_excitation Spin.Beta hole particle ki)\n",
    "              ) mos_abs\n",
    "          ) mos_b \n",
    "      ]\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 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 = particle || hole' = particle' then None else\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 && hole <> particle && hole' <> 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",
    "                  if hole = particle || hole' = particle' then None else\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 cpudet, cpuint = ref 0., ref 0. in\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",
    "        let task (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 || Determinant.degree ki kj <> exc then\n",
    "                    (0,0,0.,0.)\n",
    "                else\n",
    "                  begin\n",
    "                    let alphas = generate_alphas ki kj exc cabs lexc rexc in\n",
    "                    let det_value =\n",
    "                      let t0 = Unix.gettimeofday () in\n",
    "                      let result = compute_HaaF ki alphas kj in\n",
    "                      cpudet := !cpudet +. Unix.gettimeofday () -. t0;\n",
    "                      result\n",
    "                    in\n",
    "                    let int_value = \n",
    "                      let t0 = Unix.gettimeofday () in\n",
    "                      let result = integral_value ki kj in\n",
    "                      cpuint := !cpuint +. Unix.gettimeofday () -. t0;\n",
    "                      result\n",
    "                    in\n",
    "(*   Printf.printf \"%d %d %e %e\\n%!\" i j det_value int_value;  *)\n",
    "                    (i,j,det_value,int_value)\n",
    "                  end\n",
    "            ) det_list\n",
    "        in\n",
    "        det_list\n",
    "        |> Array.mapi (fun i ki -> (i,ki))\n",
    "        |> Array.to_list\n",
    "        |> Stream.of_list\n",
    "        |> Farm.run ~f:task\n",
    "        |> Util.stream_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 (Array.length result - 1) 0 0 0. 0. 0.\n",
    "  in\n",
    "  let error = abs_float (d -. v) in\n",
    "  if error < epsilon_float then\n",
    "  (*\n",
    "      Printf.printf \"OK: %e\\n%!\" error\n",
    "      *)\n",
    "      Printf.printf \"OK: (%d, %d) | %e %e | %e  |  cpu : %f  %f\\n%!\" i j d v error !cpudet !cpuint\n",
    "  else\n",
    "      Printf.printf \"Failed: (%d, %d) | %e %e | %e  |  cpu : %f %f\\n%!\" i j d v error !cpudet !cpuint\n",
    "\n",
    "    \n",
    "let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let array_3_init d1 d2 d3 f =\n",
    "  let result = \n",
    "      Bigarray.(Array3.create Float64 fortran_layout) d1 d2 d3\n",
    "  in\n",
    "  for k=1 to d3 do\n",
    "      for j=1 to d2 do\n",
    "          for i=1 to d1 do\n",
    "             result.{i,j,k} <- f i j k\n",
    "          done\n",
    "      done\n",
    "  done;\n",
    "  result\n",
    "\n",
    "\n",
    "let array_4_init d1 d2 d3 d4 f =\n",
    "  let result = \n",
    "      Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4 |]\n",
    "  in\n",
    "  for l=1 to d4 do\n",
    "    for k=1 to d3 do\n",
    "      for j=1 to d2 do\n",
    "        for i=1 to d1 do\n",
    "           result.{i,j,k,l} <- f i j k l\n",
    "        done\n",
    "      done\n",
    "    done\n",
    "  done;\n",
    "  result\n",
    "\n",
    "\n",
    "let array_5_init d1 d2 d3 d4 d5 f =\n",
    "  let result = \n",
    "      Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4;d5 |]\n",
    "  in\n",
    "  for m=1 to d5 do\n",
    "    for l=1 to d4 do\n",
    "      for k=1 to d3 do\n",
    "        for j=1 to d2 do\n",
    "          for i=1 to d1 do\n",
    "             result.{i,j,k,l,m} <- f i j k l m\n",
    "          done\n",
    "        done\n",
    "      done\n",
    "    done\n",
    "  done;\n",
    "  result\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": [
    "let _ = png_image \"0_1_11.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$i,j$ : orbital indices of MOs occupied in $|I\\rangle$.\n",
    "\n",
    "$\\sum_i \\sum_a \\hat{T}_{i}^{a} \\hat{T}_{a}^{i}$\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_{a} \\sum_{i}\n",
    "        \\left( h_{ia} +\n",
    "            \\sum_{j} \\langle i j || a j \\rangle +\n",
    "            \\sum_{\\bar{j}} \\langle i \\bar{j} | a \\bar{j} \\rangle \\right) \n",
    "        \\left( f_{ai} +\n",
    "            \\sum_{j}  \\left[ a j || i j \\right] + \n",
    "            \\sum_{\\bar{j}}  \\left[ a \\bar{j} | i \\bar{j} \\right] \\right) + \\\\\n",
    "\\sum_{\\bar{a}}  \\sum_{\\bar{i}}\n",
    "        \\left( h_{ia} +\n",
    "            \\sum_{j} \\langle \\bar{i} j | \\bar{a} j \\rangle +\n",
    "            \\sum_{\\bar{j}} \\langle \\bar{i} \\bar{j} || \\bar{a} \\bar{j} \\rangle \\right) \n",
    "        \\left( f_{ai} +\n",
    "            \\sum_{j}  \\left[ \\bar{a} j | \\bar{i} j \\right] + \n",
    "            \\sum_{\\bar{j}}  \\left[ \\bar{a} \\bar{j} || \\bar{i} \\bar{j} \\right] \\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",
    "   let mos_a, mos_b = mos_a ki, mos_b ki in\n",
    "   \n",
    "   sum mos_cabs (fun a ->\n",
    "      sum mos_a (fun i -> \n",
    "           (h_one i a s  +.\n",
    "               sum mos_a (fun j -> h_two i j a j s s ) +.\n",
    "               sum mos_b (fun j -> h_two i j a j s s')\n",
    "               ) *. \n",
    "           (f_one a i s  +.\n",
    "               sum mos_a (fun j -> f_two a j i j s s ) +.\n",
    "               sum mos_b (fun j -> f_two a j i j s s')\n",
    "               )\n",
    "         )\n",
    "       +.\n",
    "       sum mos_b (fun i -> \n",
    "           (h_one i a s'  +.\n",
    "               sum mos_b (fun j -> h_two i j a j s' s') +.\n",
    "               sum mos_a (fun j -> h_two i j a j s' s )\n",
    "               ) *. \n",
    "           (f_one a i s'  +.\n",
    "               sum mos_b (fun j -> f_two a j i j s' s') +.\n",
    "               sum mos_a (fun j -> f_two a j i j s' s )\n",
    "               )\n",
    "         ) \n",
    "     ) \n",
    "     \n",
    "let _ =\n",
    "  check 0 integral_value 0 1 1 1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "let m_0111_1H_1F = \n",
    "  Vec.init mo_num (fun i -> \n",
    "        sum mos_cabs (fun a ->\n",
    "            h_one i a Spin.Alfa *. f_one a i Spin.Alfa ))\n",
    "\n",
    "    \n",
    "let m_0111_1H_2Fa = \n",
    "  Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->\n",
    "    h_one i a Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Alfa  ))\n",
    "\n",
    "    \n",
    "let m_0111_1H_2Fb = \n",
    "  Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->\n",
    "    h_one i a Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Beta  ))\n",
    "\n",
    "\n",
    "let m_0111_2Ha_1F = \n",
    "  Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->\n",
    "    h_two i j a j Spin.Alfa Spin.Alfa *. f_one a i Spin.Alfa ))\n",
    "\n",
    "    \n",
    "let m_0111_2Hb_1F = \n",
    "  Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->\n",
    "    h_two i j a j Spin.Alfa Spin.Beta *. f_one a i Spin.Alfa ))\n",
    "\n",
    "\n",
    "let m_0111_2Ha_2Fa =\n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j a j Spin.Alfa Spin.Alfa *.\n",
    "      f_two a k i k Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_0111_2Hb_2Fa =\n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j a j Spin.Alfa Spin.Beta *.\n",
    "      f_two a k i k Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_0111_2Ha_2Fb =\n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j a j Spin.Alfa Spin.Alfa *.\n",
    "      f_two a k i k Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_0111_2Hb_2Fb =\n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j a j Spin.Alfa Spin.Beta *.\n",
    "      f_two a k i k Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "\n",
    "let integral_value ki kj = \n",
    "\n",
    "   let mos_a, mos_b = mos_a ki, mos_b ki in\n",
    "   \n",
    "   sum mos_a (fun i ->\n",
    "     m_0111_1H_1F.{i} +. \n",
    "     sum mos_a (fun j ->\n",
    "       m_0111_2Ha_1F.{i,j} +.\n",
    "       m_0111_1H_2Fa.{i,j} +.\n",
    "       sum mos_a (fun k -> m_0111_2Ha_2Fa.{i,j,k}) +.\n",
    "       sum mos_b (fun k -> m_0111_2Ha_2Fb.{i,j,k})\n",
    "       ) +.\n",
    "     sum mos_b (fun j ->\n",
    "       m_0111_2Hb_1F.{i,j} +.\n",
    "       m_0111_1H_2Fb.{i,j} +.\n",
    "       sum mos_a (fun k -> m_0111_2Hb_2Fa.{i,j,k}) +.\n",
    "       sum mos_b (fun k -> m_0111_2Hb_2Fb.{i,j,k})\n",
    "       )\n",
    "     )\n",
    "   +.\n",
    "   sum mos_b (fun i ->\n",
    "     m_0111_1H_1F.{i} +. \n",
    "     sum mos_b (fun j ->\n",
    "       m_0111_2Ha_1F.{i,j} +.\n",
    "       m_0111_1H_2Fa.{i,j} +.\n",
    "       sum mos_b (fun k -> m_0111_2Ha_2Fa.{i,j,k}) +.\n",
    "       sum mos_a (fun k -> m_0111_2Ha_2Fb.{i,j,k})\n",
    "       ) +.\n",
    "     sum mos_a (fun j ->\n",
    "       m_0111_2Hb_1F.{i,j} +.\n",
    "       m_0111_1H_2Fb.{i,j} +.\n",
    "       sum mos_b (fun k -> m_0111_2Hb_2Fa.{i,j,k}) +.\n",
    "       sum mos_a (fun k -> m_0111_2Hb_2Fb.{i,j,k})\n",
    "       )\n",
    "     )\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": [
    "let _ = 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",
    "$\\sum_i \\sum_a \\sum_j \\sum_k \\hat{T}_{ij}^{ak} \\hat{T}_{ak}^{ij}$\n",
    "\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",
    "   sum mos_cabs   (fun a ->\n",
    "   sum (mos_b ki) (fun j ->  \n",
    "   sum (mos_a ki) (fun i ->  \n",
    "   sum mos_virt_b (fun k ->  \n",
    "       h_two i j a k s s'  *. f_two a k i j s s'\n",
    "   )))\n",
    "   +.\n",
    "   sum (mos_a ki) (fun j ->  \n",
    "   sum (mos_b ki) (fun i ->  \n",
    "   sum mos_virt_a (fun k ->  \n",
    "       h_two i j a k s s'  *. f_two a k i j s s'\n",
    "   )))) \n",
    "   +.\n",
    "   (* Alpha-Alpha / Beta-Beta *)\n",
    "   List.fold_left (fun accu (mos_virt,mos,s) -> \n",
    "     let s' = s in accu +.\n",
    "     sum mos (fun j ->  \n",
    "     sum mos (fun i -> if i < j then 0. else\n",
    "     sum mos_virt (fun k ->  \n",
    "     sum mos_cabs (fun a -> \n",
    "       h_two i j a k s s'  *. f_two a k i j s s'\n",
    "     )))) \n",
    "   ) 0. [ (mos_virt_a,mos_a ki,Spin.Alfa) ; (mos_virt_b,mos_b ki,Spin.Beta)]\n",
    "     \n",
    "let _ =\n",
    "  check 100 integral_value 0 1 2 2\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let m_0122_Haa = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "      sum mos_cabs (fun a ->\n",
    "        h_two i j a k Spin.Alfa Spin.Alfa *. f_two a k i j Spin.Alfa Spin.Alfa\n",
    "        ) )\n",
    "      \n",
    "let m_0122_Hab = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "      sum mos_cabs (fun a ->\n",
    "         h_two i j a k Spin.Alfa Spin.Beta *. f_two a k i j Spin.Alfa Spin.Beta\n",
    "      ) )\n",
    "\n",
    "\n",
    "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",
    "   let result = \n",
    "   sum (mos_a ki) (fun i ->  \n",
    "   sum (mos_b ki) (fun j ->  \n",
    "   sum mos_virt_b (fun k ->  m_0122_Hab.{i,j,k} )))\n",
    "   +.\n",
    "   sum (mos_b ki) (fun i ->  \n",
    "   sum (mos_a ki) (fun j ->  \n",
    "   sum mos_virt_a (fun k ->  m_0122_Hab.{i,j,k} )))\n",
    "   +.\n",
    "   (* Alpha-Alpha / Beta-Beta *)\n",
    "   List.fold_left (fun accu (mos_virt,mos) -> \n",
    "     accu +.\n",
    "     sum mos (fun j ->  \n",
    "     sum mos (fun i -> if i <= j then 0. else\n",
    "     sum mos_virt (fun k -> m_0122_Haa.{i,j,k} )))\n",
    "   ) 0. [ (mos_virt_a,mos_a ki) ; (mos_virt_b,mos_b ki)]\n",
    "   in \n",
    "   result\n",
    "     \n",
    "let _ =\n",
    "  check 100 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": [
    "let _ = png_image \"0_2_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$i,j$ : orbital indices of MOs   occupied in $|I\\rangle$.\n",
    "\n",
    "$\\sum_i \\sum_a \\sum_j \\sum_b \\hat{T}_{ij}^{ab} \\hat{T}_{ab}^{ij}$\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",
    "   sum  mos_cabs  (fun a -> \n",
    "   sum  mos_cabs  (fun b ->  \n",
    "   sum (mos_a ki) (fun i ->  \n",
    "   sum (mos_b ki) (fun j ->  \n",
    "       h_two i j a b s s'  *. f_two a b i j s s'\n",
    "   ))))\n",
    "   +.\n",
    "   (* Alpha-Alpha / Beta-Beta *)\n",
    "   List.fold_left (fun accu (mos,s) -> \n",
    "     let s' = s in accu +.\n",
    "     sum mos_cabs (fun b ->  \n",
    "     sum mos_cabs (fun a -> if b > a then 0. else\n",
    "     sum   mos    (fun j ->  \n",
    "     sum   mos    (fun i -> if i < j then 0. else\n",
    "        h_two i j a b s s'  *. f_two a b i j s s'\n",
    "     ))))\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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let m_0222_Haa = \n",
    "  Mat.init_cols mo_num mo_num (fun i j ->\n",
    "      sum mos_cabs (fun a ->\n",
    "      sum mos_cabs (fun b -> if b >= a then 0. else\n",
    "        h_two i j a b Spin.Alfa Spin.Alfa *. f_two a b i j Spin.Alfa Spin.Alfa\n",
    "        ) )\n",
    "     )\n",
    "        \n",
    "let m_0222_Hab = \n",
    "  Mat.init_cols mo_num mo_num (fun i j ->\n",
    "      sum mos_cabs (fun a ->\n",
    "      sum mos_cabs (fun b ->\n",
    "        h_two i j a b Spin.Alfa Spin.Beta *. f_two a b i j Spin.Alfa Spin.Beta\n",
    "        ) )\n",
    "     )\n",
    "      \n",
    "let integral_value ki kj = \n",
    "   (* Alpha-Beta *)\n",
    "   sum (mos_a ki) (fun i ->  \n",
    "       sum (mos_b ki) (fun j ->  \n",
    "          m_0222_Hab.{i,j}  \n",
    "       ))\n",
    "   +.\n",
    "   (* Alpha-Alpha / Beta-Beta *)\n",
    "   List.fold_left (fun accu (mos) -> \n",
    "     accu +.\n",
    "     sum mos (fun j ->  \n",
    "         sum mos (fun i -> if i < j then 0. else\n",
    "           m_0222_Haa.{i,j}  \n",
    "     ))\n",
    "   ) 0. [ (mos_a ki) ; (mos_b ki)]\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": [
    "let _ = 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",
    "$\\sum_a \\hat{T}_{i}^{a} \\hat{T}_{a}^{k}$\n",
    "\n",
    "\\begin{align}\n",
    "    \\sum_{a} \n",
    "    \\left( h_{ia} +\\sum_{j} \\langle i j || a j \\rangle +\n",
    "    \\sum_{\\bar{j}} \\langle i \\bar{j} | a \\bar{j} \\rangle  \\right)\n",
    "    \\left( f_{ak} +\\sum_{j}  \\left[ a j || k j \\right] + \n",
    "    \\sum_{\\bar{j}}  \\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_i, mos_i' =\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",
    "    let mos_j, mos_j' =\n",
    "      match s with\n",
    "      | Spin.Alfa -> mos_a kj, mos_b kj\n",
    "      | Spin.Beta -> mos_b kj, mos_a kj\n",
    "    in\n",
    "\n",
    "    let result = \n",
    "       let s' = Spin.other s in\n",
    "       \n",
    "       sum mos_cabs (fun a -> \n",
    "           (h_one i a s +. \n",
    "              sum mos_i  (fun j ->  h_two i j a j s s ) +.\n",
    "              sum mos_i' (fun j ->  h_two i j a j s s') ) *.\n",
    "           (f_one a k s +. \n",
    "              sum mos_j  (fun j ->  f_two a j k j s s ) +.\n",
    "              sum mos_j' (fun j ->  f_two a j k j s s') )  \n",
    "         )\n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 0 integral_value 1 1 1 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let ki = det_I.(77)\n",
    "let kj = det_I.(83)\n",
    "\n",
    "let _ = integral_value ki kj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "let m_1111_1H_1F = \n",
    "  Mat.init_cols mo_num mo_num (fun i k -> \n",
    "        sum mos_cabs (fun a ->\n",
    "            h_one i a Spin.Alfa *. f_one a k Spin.Alfa ))\n",
    "\n",
    "\n",
    "let m_1111_1H_2Fa = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "      sum mos_cabs (fun a ->\n",
    "        h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Alfa ))\n",
    "\n",
    "    \n",
    "let m_1111_1H_2Fb = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "      sum mos_cabs (fun a ->\n",
    "        h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Beta ))\n",
    "\n",
    "\n",
    "let m_1111_2Ha_1F = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "   sum mos_cabs (fun a ->\n",
    "    h_two i j a j Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa ))\n",
    "\n",
    "    \n",
    "let m_1111_2Hb_1F = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "   sum mos_cabs (fun a ->\n",
    "    h_two i j a j Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa ))\n",
    "\n",
    "\n",
    "let m_1111_2Ha_2Fa =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j a j Spin.Alfa Spin.Alfa *.\n",
    "      f_two a l k l Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_1111_2Hb_2Fa =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j a j Spin.Alfa Spin.Beta *.\n",
    "      f_two a l k l Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_1111_2Ha_2Fb =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j a j Spin.Alfa Spin.Alfa *.\n",
    "      f_two a l k l Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_1111_2Hb_2Fb =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j a j Spin.Alfa Spin.Beta *.\n",
    "      f_two a l k l Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let integral_value_11 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_a, mos_b =\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",
    "     m_1111_1H_1F.{i,k} +. \n",
    "     sum mos_a (fun j -> if j = i then 0. else\n",
    "       m_1111_1H_2Fa.{i,j,k} +. m_1111_2Ha_1F.{i,j,k} +.\n",
    "       sum mos_a (fun l -> if l = i then 0. else m_1111_2Ha_2Fa.{i,j,k,l}) +.\n",
    "       sum mos_b (fun l -> m_1111_2Ha_2Fb.{i,j,k,l})) +.\n",
    "     sum mos_b (fun j ->\n",
    "       m_1111_1H_2Fb.{i,j,k} +. m_1111_2Hb_1F.{i,j,k} +.\n",
    "       sum mos_a (fun l -> if l = i then 0. else m_1111_2Hb_2Fa.{i,j,k,l}) +.\n",
    "       sum mos_b (fun l -> m_1111_2Hb_2Fb.{i,j,k,l}))\n",
    "    in\n",
    "    \n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "    \n",
    "let _ =\n",
    "  check 0 integral_value_11 1 1 1 1\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. 1 1 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let _ = png_image \"1_1_12.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{i}^{k}|I\\rangle$\n",
    "\n",
    "$\\sum_{j} \\sum_a \\hat{T}_{j}^{a} \\hat{T}_{ai}^{jk}$\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_{ja} +\\sum_{m} \\langle j m || a m \\rangle +\n",
    "    \\sum_{\\bar{m}} \\langle j \\bar{m} | a \\bar{m} \\rangle  \\right)\n",
    "    \\left[ a i || j k \\right] + \\\\\n",
    "    \\sum_{\\bar{a}} \\sum_{\\bar{j}} \n",
    "    \\left( h_{\\bar{j}\\bar{a}} +\\sum_{m} \\langle \\bar{j} m | \\bar{a} m \\rangle +\n",
    "    \\sum_{\\bar{m}} \\langle \\bar{j} \\bar{m} || \\bar{a} \\bar{m} \\rangle  \\right)\n",
    "    \\left[ \\bar{a} i | \\bar{j} k \\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_j, mos_j' =\n",
    "      let alfa = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      let beta  = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.beta ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.beta kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      match s with\n",
    "      | Spin.Alfa -> alfa, beta\n",
    "      | Spin.Beta -> beta, alfa\n",
    "    in\n",
    "    \n",
    "    let mos_i, mos_i' =\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",
    "       let s' = Spin.other s in\n",
    "       sum mos_cabs (fun a -> \n",
    "           sum mos_j (fun j -> \n",
    "               (h_one j a s +. \n",
    "                  sum mos_i  (fun m -> h_two j m a m s s )  +.\n",
    "                  sum mos_i' (fun m -> h_two j m a m s s')  ) *.\n",
    "               (f_two a i j k s s )\n",
    "             ) \n",
    "           +.\n",
    "           sum mos_j' (fun j -> \n",
    "               (h_one j a s +. \n",
    "                  sum mos_i  (fun m ->  h_two j m a m s' s ) +.\n",
    "                  sum mos_i' (fun m ->  h_two j m a m s' s') ) *.\n",
    "               (f_two a i j k s' s )\n",
    "             ) \n",
    "         ) \n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 100 integral_value 1 1 1 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let m_1112_1H_2Fa = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "      sum mos_cabs (fun a ->\n",
    "        h_one j a Spin.Alfa *. f_two a i j k Spin.Alfa Spin.Alfa ))\n",
    "\n",
    "    \n",
    "let m_1112_1H_2Fb = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "      sum mos_cabs (fun a ->\n",
    "        h_one j a Spin.Beta *. f_two a i j k Spin.Alfa Spin.Beta ))\n",
    "\n",
    "\n",
    "let m_1112_2Ha_2Fa =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two j m a m Spin.Alfa Spin.Alfa *.\n",
    "      f_two a i j k Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_1112_2Hb_2Fa =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two j m a m Spin.Alfa Spin.Beta *.\n",
    "      f_two a i j k Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "let m_1112_2Ha_2Fb =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two j m a m Spin.Alfa Spin.Alfa *.\n",
    "      f_two a i j k Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_1112_2Hb_2Fb =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two j m a m Spin.Alfa Spin.Beta *.\n",
    "      f_two a i j k Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let integral_value_12 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_j, mos_j' =\n",
    "      let alfa = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      let beta  = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.beta ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.beta kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      match s with\n",
    "      | Spin.Alfa -> alfa, beta\n",
    "      | Spin.Beta -> beta, alfa\n",
    "    in\n",
    "    \n",
    "    let mos_i, mos_i' =\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",
    "     sum mos_j (fun j -> \n",
    "       m_1112_1H_2Fa.{i,j,k} +. \n",
    "       sum mos_i  (fun m -> m_1112_2Ha_2Fa.{i,j,k,m}) +.\n",
    "       sum mos_i' (fun m -> m_1112_2Hb_2Fa.{i,j,k,m})) +.\n",
    "     sum mos_j' (fun j ->\n",
    "       m_1112_1H_2Fb.{i,j,k} +. \n",
    "       sum mos_i  (fun m -> m_1112_2Hb_2Fb.{i,j,k,m}) +.\n",
    "       sum mos_i' (fun m -> m_1112_2Ha_2Fb.{i,j,k,m}))\n",
    "    in\n",
    "    \n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "   \n",
    "   \n",
    "\n",
    "let _ = check 100 integral_value_12 1 1 1 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. 1 1 21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "let _ = png_image \"1_1_21.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{i}^{k}|I\\rangle$\n",
    "\n",
    "$\\sum_{j} \\sum_a \\hat{T}_{ij}^{ka} \\hat{T}_{a}^{j}$\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",
    "    \\langle i j || k a \\rangle \n",
    "    \\left( f_{aj} +\\sum_{m} \\left[ m a || m j \\right] +\n",
    "    \\sum_{\\bar{m}} \\left[ \\bar{m} a | \\bar{m} j \\right] \\right) + \\\\\n",
    "    \\sum_{\\bar{a}} \\sum_{\\bar{j}} \n",
    "    \\langle i \\bar{j} | k \\bar{a} \\rangle\n",
    "    \\left( f_{\\bar{a}\\bar{j}} +\\sum_{m} \\left[ m \\bar{a} | m \\bar{j} \\right] +\n",
    "    \\sum_{\\bar{m}} \\left[ \\bar{m} \\bar{a} || \\bar{m} \\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_j, mos_j' =\n",
    "      let alfa = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      let beta  = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.beta ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.beta kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      match s with\n",
    "      | Spin.Alfa -> alfa, beta\n",
    "      | Spin.Beta -> beta, alfa\n",
    "    in\n",
    "    \n",
    "    let mos_i, mos_i' =\n",
    "      match s with\n",
    "      | Spin.Alfa -> mos_a kj, mos_b kj\n",
    "      | Spin.Beta -> mos_b kj, mos_a kj\n",
    "    in\n",
    "\n",
    "    let result = \n",
    "       let s' = Spin.other s in\n",
    "       sum mos_cabs (fun a -> \n",
    "       sum mos_j    (fun j -> \n",
    "               (h_two i j k a s s ) *.\n",
    "               (f_one a j s  +. \n",
    "                  sum mos_i  (fun m -> f_two m a m j s  s) +.\n",
    "                  sum mos_i' (fun m -> f_two m a m j s' s) ) \n",
    "             ) \n",
    "       +.\n",
    "       sum mos_j' (fun j -> \n",
    "               (h_two i j k a s s') *.\n",
    "               (f_one a j s' +. \n",
    "                  sum mos_i  (fun m -> f_two m a m j s  s') +.\n",
    "                  sum mos_i' (fun m -> f_two m a m j s' s') ) \n",
    "       )) \n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 100 integral_value 1 1 2 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let m_1121_2Ha_1F = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "      sum mos_cabs (fun a ->\n",
    "        h_two i j k a Spin.Alfa Spin.Alfa *. f_one a j Spin.Alfa))\n",
    "\n",
    "    \n",
    "let m_1121_2Hb_1F = \n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "      sum mos_cabs (fun a ->\n",
    "        h_two i j k a Spin.Alfa Spin.Beta *. f_one a j Spin.Beta))\n",
    "\n",
    "\n",
    "let m_1121_2Ha_2Fa =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j k a Spin.Alfa Spin.Alfa *.\n",
    "      f_two m a m j Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_1121_2Hb_2Fa =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j k a Spin.Alfa Spin.Beta *.\n",
    "      f_two m a m j Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "let m_1121_2Ha_2Fb =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j k a Spin.Alfa Spin.Alfa *.\n",
    "      f_two m a m j Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let m_1121_2Hb_2Fb =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j k a Spin.Alfa Spin.Beta *.\n",
    "      f_two m a m j Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let integral_value_21 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_j, mos_j' =\n",
    "      let alfa = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      let beta  = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.beta ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.beta kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      match s with\n",
    "      | Spin.Alfa -> alfa, beta\n",
    "      | Spin.Beta -> beta, alfa\n",
    "    in\n",
    "    \n",
    "    let mos_i, mos_i' =\n",
    "      match s with\n",
    "      | Spin.Alfa -> mos_a kj, mos_b kj\n",
    "      | Spin.Beta -> mos_b kj, mos_a kj\n",
    "    in\n",
    "   \n",
    "    let result = \n",
    "     sum mos_j (fun j -> \n",
    "       m_1121_2Ha_1F.{i,j,k} +. \n",
    "       sum mos_i  (fun m -> m_1121_2Ha_2Fa.{i,j,k,m}) +.\n",
    "       sum mos_i' (fun m -> m_1121_2Ha_2Fb.{i,j,k,m})) +.\n",
    "     sum mos_j' (fun j ->\n",
    "       m_1121_2Hb_1F.{i,j,k} +. \n",
    "       sum mos_i  (fun m -> m_1121_2Hb_2Fb.{i,j,k,m}) +.\n",
    "       sum mos_i' (fun m -> m_1121_2Hb_2Fa.{i,j,k,m}))\n",
    "    in\n",
    "    \n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "   \n",
    "   \n",
    "\n",
    "let _ = check 100 integral_value_21 1 1 2 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. 1 1 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let _ = png_image \"1_1_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{i}^{k}|I\\rangle$\n",
    "\n",
    "$\\sum_{j} \\sum_{a} \\sum_{m}\n",
    "    \\hat{T}_{ij}^{ma} \\hat{T}_{ma}^{kj} - \n",
    "    \\hat{T}_{lj}^{ka} \\hat{T}_{ia}^{lj} + \\\\\n",
    "\\sum_{\\bar{j}} \\sum_{\\bar{a}} \\sum_{m}\n",
    "    \\hat{T}_{i\\bar{j}}^{m\\bar{a}} \\hat{T}_{m\\bar{a}}^{k\\bar{j}} - \n",
    "    \\hat{T}_{l\\bar{j}}^{k\\bar{a}} \\hat{T}_{i\\bar{a}}^{l\\bar{j}} +\n",
    "\\sum_{\\bar{j}} \\sum_{a} \\sum_{\\bar{m}}\n",
    "    \\hat{T}_{i\\bar{j}}^{a\\bar{m}} \\hat{T}_{a\\bar{m}}^{k\\bar{j}} \n",
    "$\n",
    "\n",
    "$j,l$ : orbital indices of MOs occupied both in $|I\\rangle$  and $|J\\rangle$.\n",
    "\n",
    "$m$ : orbital indices of virtual MOs, unoccupied both in $|I\\rangle$  and $|J\\rangle$.\n",
    "\n",
    "$a$ : orbital indices of CABS MOs\n",
    "\n",
    "\\begin{align}\n",
    "    \\sum_{a} \\sum_{j} \\sum_{m}\n",
    "      \\langle i j || m a \\rangle  \\left[ m a || k j \\right]  - \\\\\n",
    "    \\sum_{a} \\sum_{j} \\sum_{l}\n",
    "      \\langle l j || k a \\rangle  \\left[ i a || l j \\right]  + \\\\\n",
    "    \\sum_{\\bar{a}} \\sum_{\\bar{j}} \\sum_{m}\n",
    "      \\langle i \\bar{j} | m \\bar{a} \\rangle  \\left[ m \\bar{a} | k \\bar{j} \\right]  + \\\\\n",
    "    \\sum_{\\bar{a}} \\sum_{\\bar{j}} \\sum_{l}\n",
    "      \\langle l \\bar{j} | k \\bar{a} \\rangle  \\left[ i \\bar{a} | l \\bar{j} \\right] - \\\\\n",
    "    \\sum_{\\bar{a}} \\sum_{\\bar{j}} \\sum_{\\bar{m}}\n",
    "      \\langle i \\bar{j} | a \\bar{m} \\rangle  \\left[ a \\bar{m} | k \\bar{j} \\right]  \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(*\n",
    "let ki = det_I.(6)\n",
    "let kj = det_I.(58)\n",
    "\n",
    "let _ = \n",
    "    Format.printf \"|I> -> |J> : %a |\\n@.\" Excitation.pp (Excitation.of_det ki kj) ;\n",
    "generate_alphas ki kj 1 1 2 2\n",
    "|> Array.of_list\n",
    "|> Array.mapi (fun kk alpha -> \n",
    "\n",
    "(*\n",
    "let _ = Determinant.to_lists ki\n",
    "let _ = Determinant.to_lists alpha\n",
    "let _ = Determinant.to_lists kj\n",
    "\n",
    "   \n",
    "let _ = \n",
    "    Format.printf \"|I> -> |J> : %a\\n@.\" Excitation.pp (Excitation.of_det ki kj);\n",
    "    Format.printf \"|I> -> |a> : %a\\n@.\" Excitation.pp (Excitation.of_det ki alpha);\n",
    "    Format.printf \"|a> -> |J> : %a\\n@.\" Excitation.pp (Excitation.of_det alpha kj)\n",
    "*)\n",
    "\n",
    "    let integral_value ki kj = \n",
    "    \n",
    "      let exc0 = Array.init (aux_num+1) (fun _ -> [|\"-\";\"-\"|]) in\n",
    "      let exc1 = Array.init (aux_num+1) (fun _ -> [|\"-\";\"-\"|]) in \n",
    "      let exc2 = Array.init (aux_num+1) (fun _ -> [|\"-\";\"-\"|]) in \n",
    " \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",
    "      let spin = function\n",
    "        | Spin.Alfa -> 0\n",
    "        | _ -> 1\n",
    "      in\n",
    "        exc0.(i).(spin s ) <- \"i\" \n",
    "      ; exc0.(k).(spin s ) <- \"k\" \n",
    "      ;\n",
    "      let s0 = s in\n",
    "    \n",
    "      let i, j, k, l, s, s', p1 =\n",
    "          match Excitation.of_det ki alpha with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      \n",
    "      if exc0.(i).(spin s ) =  \"-\" then exc0.(i).(spin s ) <- \"n\";\n",
    "      if exc0.(j).(spin s') =  \"-\" then exc0.(j).(spin s') <- \"n\";\n",
    "      if exc0.(k).(spin s ) =  \"-\" then exc0.(k).(spin s ) <- if k > mo_num then \"a\" else \"m\";\n",
    "      if exc0.(l).(spin s') =  \"-\" then exc0.(l).(spin s') <- if l > mo_num then \"a\" else \"m\";\n",
    "      \n",
    "      let string_h = \n",
    "          Printf.sprintf \"h_two %s %s %s %s %s %s *.  \" \n",
    "            exc0.(i).(spin s )\n",
    "            exc0.(j).(spin s')\n",
    "            exc0.(k).(spin s )\n",
    "            exc0.(l).(spin s')\n",
    "            (if s  = s0 then \"s \" else \"s'\")\n",
    "            (if s' = s0 then \"s \" else \"s'\")\n",
    "            (*\n",
    "            (if exc0.(i).(spin s ) = \"i\" || exc0.(k).(spin s ) = \"k\" then \"s  \" else \n",
    "             if exc0.(i).(spin s ) = \"j\" || exc0.(k).(spin s ) = \"l\" then \"s' \" else \"s''\")\n",
    "            (if exc0.(j).(spin s') = \"i\" || exc0.(l).(spin s') = \"k\" then \"s  \" else \n",
    "             if exc0.(j).(spin s') = \"j\" || exc0.(l).(spin s') = \"l\" then \"s' \" else \"s''\")\n",
    "             *)\n",
    "      in \n",
    "      \n",
    "      let i, j, k, l, s, s', p2 =\n",
    "          match Excitation.of_det alpha kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      \n",
    "      let string_f = \n",
    "          Printf.sprintf \"f_two %s %s %s %s %s %s\" \n",
    "            exc0.(i).(spin s )\n",
    "            exc0.(j).(spin s')\n",
    "            exc0.(k).(spin s )\n",
    "            exc0.(l).(spin s')\n",
    "            (*\n",
    "            (if exc0.(i).(spin s ) = \"i\" || exc0.(k).(spin s ) = \"k\" then \"s  \" else \n",
    "             if exc0.(i).(spin s ) = \"j\" || exc0.(k).(spin s ) = \"l\" then \"s' \" else \n",
    "             if exc0.(i).(spin s ) = \"n\" || exc0.(k).(spin s ) = \"n\" then \"s''\" else\n",
    "             if s  = s0 then \"s\" else \"s'\")\n",
    "            (if exc0.(j).(spin s') = \"i\" || exc0.(l).(spin s') = \"k\" then \"s  \" else \n",
    "             if exc0.(j).(spin s') = \"j\" || exc0.(l).(spin s') = \"l\" then \"s' \" else \n",
    "             if exc0.(j).(spin s') = \"n\" || exc0.(l).(spin s') = \"n\" then \"s''\" else\n",
    "             if s' = s0 then \"s\" else \"s'\")\n",
    "             *)\n",
    "            (if s  = s0 then \"s \" else \"s'\")\n",
    "            (if s' = s0 then \"s \" else \"s'\")\n",
    "      in\n",
    "      \n",
    "    Format.printf \"|I> -> |a> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det ki alpha) string_h ;\n",
    "    Format.printf \"|a> -> |J> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det alpha kj) string_f ;\n",
    "    \n",
    "    let sign = \n",
    "        if Phase.add p1 p2 = phase then \"+. \" else \"-. \"\n",
    "    in\n",
    "sign ^ string_h ^ string_f\n",
    "(*\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",
    "           let k=4 and l=5 in\n",
    "           let i=1 and j=2 in\n",
    "           let n=2 and a=7 in\n",
    "           let s   = Spin.Alfa\n",
    "           and s'  = Spin.Alfa\n",
    "           and s'' = Spin.Beta\n",
    "           in\n",
    "           (* \n",
    "           h_two j n l a s' s'' *. f_two a i n k s s''\n",
    "           *)\n",
    "           h_two i n k a s s'' *. f_two j a l n s s''\n",
    "        in\n",
    "\n",
    "        match phase with\n",
    "        | Phase.Pos -> result\n",
    "        | Phase.Neg -> -. result\n",
    "        *)\n",
    "\n",
    "in\n",
    "integral_value ki kj\n",
    "(*\n",
    "let a = (compute_HaaF ki [alpha] kj)\n",
    "and b = (integral_value ki kj)\n",
    "in\n",
    "if kk = 31 then\n",
    "    Format.printf \"%6d %e %e@.@.\" (kk) a b\n",
    "    *)\n",
    ")\n",
    "|> Array.to_list\n",
    "|> List.sort_uniq compare\n",
    "\n",
    "let _ = compute_HaaF ki (generate_alphas ki kj 1 1 2 2) kj\n",
    "\n",
    "let _ = integral_value ki kj\n",
    "\n",
    "*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "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",
    "   (* 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",
    "   let mos_virt, mos_virt' =\n",
    "      match s with\n",
    "      | Spin.Alfa -> mos_virt_a, mos_virt_b\n",
    "      | Spin.Beta -> mos_virt_b, mos_virt_a\n",
    "   in\n",
    "    \n",
    "   let mos_j, mos_j' =\n",
    "      let alfa = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      let beta  = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.beta ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.beta kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      match s with\n",
    "      | Spin.Alfa -> alfa, beta\n",
    "      | Spin.Beta -> beta, alfa\n",
    "   in\n",
    "   let result = \n",
    "       let s' = Spin.other s in\n",
    "        sum mos_cabs (fun a -> \n",
    "        sum mos_j    (fun n -> \n",
    "        sum mos_virt (fun m -> \n",
    "             h_two n i a m s  s  *.  f_two m a k n s  s \n",
    "        ))) +.\n",
    "        sum mos_cabs (fun a' -> \n",
    "        sum mos_j'   (fun n' -> \n",
    "        sum mos_virt (fun m -> \n",
    "             h_two i n' m a' s  s' *.  f_two m a' k n' s  s'\n",
    "        ))) +.\n",
    "        sum mos_cabs  (fun a -> \n",
    "        sum mos_j'    (fun n' -> \n",
    "        sum mos_virt' (fun m' -> \n",
    "             h_two i n' a m' s  s' *.  f_two a m' k n' s  s'\n",
    "        ))) -. \n",
    "        sum mos_cabs (fun a -> \n",
    "        sum mos_j    (fun n -> \n",
    "        sum mos_j    (fun l -> if l <=n then 0. else\n",
    "             h_two n l a k s  s  *.  f_two i a l n s  s \n",
    "        ))) -.  \n",
    "        sum mos_cabs (fun a' -> \n",
    "        sum mos_j'   (fun n' -> \n",
    "        sum mos_j    (fun l -> \n",
    "             h_two l n' k a' s s'  *.  f_two i a' l n' s s'\n",
    "        ))) \n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "let _ = check 100 integral_value 1 1 2 2\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "let m_1122_va =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two j i a m Spin.Alfa Spin.Alfa *.\n",
    "      f_two m a k j Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "let m_1122_v2 =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j m a Spin.Alfa Spin.Beta *.\n",
    "      f_two m a k j Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "let m_1122_v3 =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two i j a m Spin.Alfa Spin.Beta *.\n",
    "      f_two a m k j Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "let m_1122_oa =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two j l a k Spin.Alfa Spin.Alfa *.\n",
    "      f_two i a l j Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "  )\n",
    "\n",
    "let m_1122_o =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "      h_two l j k a Spin.Alfa Spin.Beta *.\n",
    "      f_two i a l j Spin.Alfa Spin.Beta \n",
    "    )\n",
    "  )\n",
    "\n",
    "\n",
    "let integral_value_22 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",
    "   (* 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",
    "   let mos_virt, mos_virt' =\n",
    "      match s with\n",
    "      | Spin.Alfa -> mos_virt_a, mos_virt_b\n",
    "      | Spin.Beta -> mos_virt_b, mos_virt_a\n",
    "   in\n",
    "    \n",
    "   let mos_j, mos_j' =\n",
    "      let alfa = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      let beta  = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.beta ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.beta kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      match s with\n",
    "      | Spin.Alfa -> alfa, beta\n",
    "      | Spin.Beta -> beta, alfa\n",
    "   in\n",
    "   \n",
    "    let result = \n",
    "        sum mos_j    (fun j -> \n",
    "        sum mos_virt (fun m -> \n",
    "          m_1122_va.{i,j,k,m}\n",
    "        )) +.\n",
    "        sum mos_j'   (fun j -> \n",
    "        sum mos_virt (fun m -> \n",
    "          m_1122_v2.{i,j,k,m}\n",
    "        ))   +.\n",
    "        \n",
    "        sum mos_j'    (fun j -> \n",
    "        sum mos_virt' (fun m -> \n",
    "          m_1122_v3.{i,j,k,m}\n",
    "        ))   -.\n",
    "        sum mos_j    (fun j -> \n",
    "        sum mos_j    (fun l -> if l <= j then 0. else\n",
    "          m_1122_oa.{i,j,k,l}\n",
    "        )) -. \n",
    "        sum mos_j'   (fun j -> \n",
    "        sum mos_j    (fun l -> \n",
    "          m_1122_o.{i,j,k,l}\n",
    "        )) \n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 100 integral_value_22 1 1 2 2\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let ki = det_I.(77) and kj = det_I.(83)\n",
    "\n",
    "let _ =\n",
    "integral_value_11 ki kj\n",
    "+. integral_value_12 ki kj\n",
    "+. integral_value_21 ki kj\n",
    "+. integral_value_22 ki kj "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. 1 2 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let _ = 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",
    "$\\sum_{j} \\sum_{a} \\sum_{b}\n",
    "    \\hat{T}_{ij}^{ab} \\hat{T}_{ab}^{kj} +\n",
    "\\sum_{\\bar{j}} \\sum_{a} \\sum_{\\bar{b}}\n",
    "    \\hat{T}_{i\\bar{j}}^{a\\bar{b}} \\hat{T}_{a\\bar{b}}^{k\\bar{j}} \n",
    "$\n",
    "\n",
    "\n",
    "$$\n",
    "\\sum_{j} \\sum_{b}\\sum_{a<b}\n",
    "    \\langle i j || a b \\rangle \\left[ ab || kj \\right] + \n",
    "\\sum_{\\bar{j}} \\sum_{\\bar{b}}\\sum_{a}\n",
    "    \\langle i \\bar{j} | a \\bar{b} \\rangle \\left[ a \\bar{b} | k \\bar{j} \\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",
    "   \n",
    "\n",
    "    let result = \n",
    "       (* Alpha-Beta *)\n",
    "       let s' = Spin.other s in\n",
    "       sum mos_cabs (fun b -> \n",
    "       sum mos_cabs (fun a -> \n",
    "       sum mos'     (fun j -> h_two h j a b s s'  *. f_two a b p j s s')\n",
    "       ))\n",
    "       +.\n",
    "       (* Alpha-Alpha / Beta-Beta *)\n",
    "       sum mos_cabs (fun b -> \n",
    "       sum mos_cabs (fun a -> if b >= a then 0. else\n",
    "       sum mos  (fun j -> h_two h j a b s s  *. f_two a b p j s s)\n",
    "       )) \n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 100 integral_value 1 2 2 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "let m_1222a =\n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "    sum mos_cabs (fun a ->\n",
    "    sum mos_cabs (fun b -> if b > a then 0. else\n",
    "      h_two i j a b Spin.Alfa Spin.Alfa *.\n",
    "      f_two a b k j Spin.Alfa Spin.Alfa \n",
    "    )\n",
    "    )\n",
    "  )\n",
    "\n",
    "let m_1222 =\n",
    "  array_3_init mo_num mo_num mo_num (fun i j k ->\n",
    "    sum mos_cabs (fun a ->\n",
    "    sum mos_cabs (fun b ->\n",
    "      h_two i j a b Spin.Alfa Spin.Beta *.\n",
    "      f_two a b k j Spin.Alfa Spin.Beta \n",
    "    )\n",
    "    )\n",
    "  )\n",
    "  \n",
    "\n",
    "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",
    "   (* MOs unoccupied in both I and J *)\n",
    "    let mos_j, mos_j' =\n",
    "      let alfa = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      let beta  = \n",
    "        let i = Spindeterminant.bitstring @@ Determinant.beta ki in\n",
    "        let j = Spindeterminant.bitstring @@ Determinant.beta kj in\n",
    "        Bitstring.to_list (Bitstring.logand i j)\n",
    "      in\n",
    "      match s with\n",
    "      | Spin.Alfa -> alfa, beta\n",
    "      | Spin.Beta -> beta, alfa\n",
    "    in\n",
    "\n",
    "    let result = \n",
    "        sum mos_j    (fun j ->  m_1222a.{i,j,k} ) +.\n",
    "        sum mos_j'   (fun j ->  m_1222.{i,j,k} )\n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 100 integral_value 1 2 2 2\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9. 2 1 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let _ = png_image \"2_1_12.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
    "\n",
    "$m$ : orbital indices of MOs unoccupied both in $|I\\rangle$  and $|J\\rangle$.\n",
    "\n",
    "$n$ : orbital indices of MOs occupied in $|I\\rangle$\n",
    "\n",
    "$\\sum_{a} \n",
    "    \\hat{T}_{i}^{a} \\hat{T}_{aj}^{kl} +\n",
    "    \\hat{T}_{j}^{a} \\hat{T}_{ia}^{kl} \n",
    "$\n",
    "\n",
    "\n",
    "$$\n",
    "\\sum_{a}\n",
    "    \\left( h_{ia} + \\sum_{n} \\langle i n || a n \\rangle + \\sum_{\\bar{n}} \\langle i \\bar{n} | a \\bar{n} \\rangle \\right) \\left[ a j || k l \\right] + \\\\\n",
    "\\sum_{a}\n",
    "    \\left( h_{ja} + \\sum_{n} \\langle j n || a n \\rangle + \\sum_{\\bar{n}} \\langle j \\bar{n} | a \\bar{n} \\rangle \\right) \\left[ i a || k l \\right] \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "    let i, j, k, l, s, s', phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, 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 mos2, mos2' =\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",
    "       sum mos_cabs (fun a -> \n",
    "           let s'' = Spin.other s in\n",
    "             ( h_one i a s +.\n",
    "               sum mos  (fun n -> h_two i n a n s s  ) +.\n",
    "               sum mos' (fun n -> h_two i n a n s s'') \n",
    "             ) *. f_two a j k l s s'\n",
    "             ) +.\n",
    "       sum mos_cabs (fun a -> \n",
    "           let s'' = Spin.other s' in\n",
    "             ( h_one j a s' +.\n",
    "               sum mos2 (fun n -> h_two j n a n s' s' ) +.\n",
    "               sum mos2'  (fun n -> h_two j n a n s' s'') \n",
    "             ) *. f_two i a k l s s'\n",
    "             ) \n",
    "        in\n",
    "\n",
    "        match phase with\n",
    "        | Phase.Pos -> result\n",
    "        | Phase.Neg -> -. result\n",
    "\n",
    "\n",
    "let _ = check 100 integral_value 2 1 1 2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
    "\n",
    "$m$ : orbital indices of MOs unoccupied both in $|I\\rangle$  and $|J\\rangle$.\n",
    "\n",
    "$n$ : orbital indices of MOs occupied in $|I\\rangle$\n",
    "\n",
    "$\\sum_{a} \n",
    "    \\hat{T}_{i}^{a} \\hat{T}_{aj}^{kl} +\n",
    "    \\hat{T}_{j}^{a} \\hat{T}_{ia}^{kl} \n",
    "$\n",
    "\n",
    "\n",
    "$$\n",
    "\\sum_{a}\n",
    "    \\left( h_{ia} + \\sum_{n} \\langle i n || a n \\rangle + \\sum_{\\bar{n}} \\langle i \\bar{n} | a \\bar{n} \\rangle \\right) \\left[ a j || k l \\right] + \\\\\n",
    "\\sum_{a}\n",
    "    \\left( h_{ja} + \\sum_{n} \\langle j n || a n \\rangle + \\sum_{\\bar{n}} \\langle j \\bar{n} | a \\bar{n} \\rangle \\right) \\left[ i a || k l \\right] \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "let m_2112_1H_2Fa =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_one i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +.\n",
    "        h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Alfa )\n",
    "    )\n",
    "\n",
    "let m_2112_1H_2Fb =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_one i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.\n",
    "        h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Beta)\n",
    "    )\n",
    "\n",
    "let m_2112_2Ha_2Fa =\n",
    "  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i n a n Spin.Alfa Spin.Alfa *.\n",
    "        f_two a j k l Spin.Alfa Spin.Alfa )\n",
    "    )\n",
    "    \n",
    "let m_2112_2Hb_2Fa =\n",
    "  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i n a n Spin.Alfa Spin.Beta *.\n",
    "        f_two a j k l Spin.Alfa Spin.Alfa )\n",
    "    )\n",
    "    \n",
    "let m_2112_2Ha_2Fb =\n",
    "  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i n a n Spin.Alfa Spin.Alfa *.\n",
    "        f_two a j k l Spin.Alfa Spin.Beta )\n",
    "    )\n",
    "    \n",
    "let m_2112_2Hb_2Fb =\n",
    "  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i n a n Spin.Alfa Spin.Beta *.\n",
    "        f_two a j k l Spin.Alfa Spin.Beta )\n",
    "    )\n",
    "    \n",
    "let integral_value ki kj = \n",
    "    let i, j, k, l, s, s', phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "    in\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",
    "       if s = s' then\n",
    "         m_2112_1H_2Fa.{i,j,k,l} +. \n",
    "         sum mos  (fun n ->\n",
    "             m_2112_2Ha_2Fa.{i,j,k,l,n} +. m_2112_2Ha_2Fa.{j,i,l,k,n}\n",
    "         ) +.\n",
    "         sum mos' (fun n ->\n",
    "             m_2112_2Hb_2Fa.{i,j,k,l,n} +. m_2112_2Hb_2Fa.{j,i,l,k,n}\n",
    "         ) \n",
    "       else\n",
    "         m_2112_1H_2Fb.{i,j,k,l} +.\n",
    "         sum mos  (fun n ->\n",
    "             m_2112_2Ha_2Fb.{i,j,k,l,n} +. m_2112_2Hb_2Fb.{j,i,l,k,n}\n",
    "         ) +.\n",
    "         sum mos' (fun n ->\n",
    "             m_2112_2Hb_2Fb.{i,j,k,l,n} +. m_2112_2Ha_2Fb.{j,i,l,k,n}\n",
    "         ) \n",
    "     in\n",
    "\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "let _ = check 100 integral_value 2 1 1 2\n",
    "\n",
    "let ki = det_I.(2)\n",
    "let kj = det_I.(33)\n",
    "\n",
    "let _ = integral_value ki kj\n",
    "let _ = \n",
    "    let alphas =  generate_alphas ki kj 2 1 1 2\n",
    "    in compute_HaaF ki alphas kj\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10. 2 1 21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let _ = png_image \"2_1_21.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
    "\n",
    "$m$ : orbital indices of MOs unoccupied both in $|I\\rangle$  and $|J\\rangle$.\n",
    "\n",
    "$n$ : orbital indices of MOs occupied in $|J\\rangle$\n",
    "\n",
    "$\\sum_{a} \n",
    "    \\hat{T}_{ij}^{al} \\hat{T}_{a}^{k} +\n",
    "    \\hat{T}_{ij}^{ka} \\hat{T}_{a}^{l} +\n",
    "    \\hat{T}_{i\\bar{j}}^{a\\bar{l}} \\hat{T}_{a}^{k} +\n",
    "    \\hat{T}_{i\\bar{j}}^{k\\bar{a}} \\hat{T}_{\\bar{a}}^{\\bar{l}}\n",
    "$\n",
    "\n",
    "\n",
    "$$\n",
    "\\sum_{a}\n",
    "    \\langle i j || a l \\rangle \\left( f_{ak} + \\sum_{n} \\left[ a n || k n \\right] + \\sum_{\\bar{n}} \\left[ a \\bar{n} | k \\bar{n} \\right] \\right) + \\\\\n",
    "\\sum_{a}\n",
    "    \\langle i j || k a \\rangle \\left( f_{al} + \\sum_{n} \\left[ n a || n l \\right] + \\sum_{\\bar{n}} \\left[ \\bar{n} a | \\bar{n} l \\right] \\right) \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "    let i, j, k, l, s, s', phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "    in\n",
    "\n",
    "        let result = \n",
    "        \n",
    "        let mos, mos', s'' =\n",
    "          match s with\n",
    "          | Spin.Alfa -> mos_a kj, mos_b kj, Spin.Beta\n",
    "          | Spin.Beta -> mos_b kj, mos_a kj, Spin.Alfa\n",
    "        in\n",
    "        sum mos_cabs (fun a -> \n",
    "             h_two i j a l s s' *.\n",
    "             ( f_one a k s +.\n",
    "               sum mos  (fun n -> f_two a n k n s s) +.\n",
    "               sum mos' (fun n -> f_two a n k n s s'') \n",
    "             ) ) +.\n",
    "\n",
    "        let mos, mos', s'' =\n",
    "          match s' with\n",
    "          | Spin.Alfa -> mos_a kj, mos_b kj, Spin.Beta\n",
    "          | Spin.Beta -> mos_b kj, mos_a kj, Spin.Alfa\n",
    "        in\n",
    "        sum mos_cabs (fun a -> \n",
    "             h_two j i a k s' s *.\n",
    "             ( f_one a l s' +.\n",
    "               sum mos  (fun n -> f_two a n l n s' s') +.\n",
    "               sum mos' (fun n -> f_two a n l n s' s'') \n",
    "             ) \n",
    "             ) \n",
    "        in\n",
    "\n",
    "        match phase with\n",
    "        | Phase.Pos -> result\n",
    "        | Phase.Neg -> -. result\n",
    "\n",
    "\n",
    "let _ = check 100 integral_value 2 1 2 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "let m_2121_2Ha_1F =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i j a l Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.\n",
    "        h_two i j k a Spin.Alfa Spin.Alfa *. f_one a l Spin.Alfa)\n",
    "    )\n",
    "\n",
    "let m_2121_2Hb_1F =\n",
    "  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i j a l Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.\n",
    "        h_two i j k a Spin.Alfa Spin.Beta *. f_one a l Spin.Alfa)\n",
    "    )\n",
    "\n",
    "let m_2121_2Ha_2Fa =\n",
    "  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i j a l Spin.Alfa Spin.Alfa *.\n",
    "        f_two a n k n Spin.Alfa Spin.Alfa )\n",
    "    )\n",
    "    \n",
    "let m_2121_2Hb_2Fa =\n",
    "  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i j a l Spin.Alfa Spin.Beta *.\n",
    "        f_two a n k n Spin.Alfa Spin.Alfa )\n",
    "    )\n",
    "    \n",
    "let m_2121_2Ha_2Fb =\n",
    "  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i j a l Spin.Alfa Spin.Alfa *.\n",
    "        f_two a n k n Spin.Alfa Spin.Beta )\n",
    "    )\n",
    "    \n",
    "let m_2121_2Hb_2Fb =\n",
    "  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->\n",
    "    sum mos_cabs (fun a ->\n",
    "        h_two i j a l Spin.Alfa Spin.Beta *.\n",
    "        f_two a n k n Spin.Alfa Spin.Beta )\n",
    "    )\n",
    "    \n",
    "let integral_value ki kj = \n",
    "    let i, j, k, l, s, s', phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "    in\n",
    "    let mos, mos' =\n",
    "       match s with\n",
    "       | Spin.Alfa -> mos_a kj, mos_b kj\n",
    "       | Spin.Beta -> mos_b kj, mos_a kj\n",
    "     in\n",
    "     \n",
    "     let result = \n",
    "       if s = s' then\n",
    "         m_2121_2Ha_1F.{i,j,k,l} +. \n",
    "         sum mos  (fun n ->\n",
    "             m_2121_2Ha_2Fa.{i,j,k,l,n} +. m_2121_2Ha_2Fa.{j,i,l,k,n}\n",
    "         ) +.\n",
    "         sum mos' (fun n ->\n",
    "             m_2121_2Ha_2Fb.{i,j,k,l,n} +. m_2121_2Ha_2Fb.{j,i,l,k,n}\n",
    "         ) \n",
    "       else\n",
    "         m_2121_2Hb_1F.{i,j,k,l} +.\n",
    "         sum mos  (fun n ->\n",
    "             m_2121_2Hb_2Fa.{i,j,k,l,n} +. m_2121_2Hb_2Fb.{j,i,l,k,n}\n",
    "         ) +.\n",
    "         sum mos' (fun n ->\n",
    "             m_2121_2Hb_2Fb.{i,j,k,l,n} +. m_2121_2Hb_2Fa.{j,i,l,k,n}\n",
    "         ) \n",
    "     in\n",
    "\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 100 integral_value 2 1 2 1\n",
    "\n",
    "(*\n",
    "let ki = det_I.(2)\n",
    "let kj = det_I.(33)\n",
    "\n",
    "let _ = integral_value ki kj\n",
    "let _ = \n",
    "    let alphas =  generate_alphas ki kj 2 1 2 1\n",
    "    in compute_HaaF ki alphas kj\n",
    "*)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 11. 2 1 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let _ = png_image \"2_1_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
    "\n",
    "$m$ : orbital indices of MOs unoccupied both in $|I\\rangle$  and $|J\\rangle$.\n",
    "\n",
    "$n$ : orbital indices of MOs occupied in both $|I\\rangle$ and $|J\\rangle$ .\n",
    "\n",
    "$\\sum_{a} \n",
    "    \\hat{T}_{ij}^{am} \\hat{T}_{am}^{kl} +\n",
    "    \\hat{T}_{ij}^{ma} \\hat{T}_{ma}^{kl} -\n",
    "    \\hat{T}_{in}^{al} \\hat{T}_{aj}^{kn} -\n",
    "    \\hat{T}_{nj}^{ka} \\hat{T}_{ia}^{nl} +\n",
    "    \\sum_{n}\n",
    "    \\hat{T}_{i\\bar{n}}^{k\\bar{a}} \\hat{T}_{\\bar{a}j}^{\\bar{n}l}\n",
    "$\n",
    "\n",
    "\n",
    "$$\n",
    "\\sum_{a}\n",
    "\\sum_{m} \\langle i j || a m \\rangle \\left[ a m || k l \\right] +  \\langle i j || m a \\rangle \\left[ m a || k l \\right] - \\\\\n",
    "\\sum_{n} \\langle i n || a l \\rangle \\left[ a j || k n \\right] +  \\langle n j || k a \\rangle \\left[ i a || n l \\right] \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "    let i, j, k, l, s, s', phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "    in\n",
    "\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",
    "    let result = \n",
    "        \n",
    "        let mos_virt, mos_virt' =\n",
    "          match s with\n",
    "          | Spin.Alfa -> mos_virt_a, mos_virt_b\n",
    "          | Spin.Beta -> mos_virt_b, mos_virt_a\n",
    "        in\n",
    "        \n",
    "        let mos, mos' =\n",
    "          let alfa = \n",
    "            let i = Spindeterminant.bitstring @@ Determinant.alfa ki in\n",
    "            let j = Spindeterminant.bitstring @@ Determinant.alfa kj in\n",
    "            Bitstring.to_list (Bitstring.logand i j)\n",
    "          in\n",
    "          let beta  =  \n",
    "            let i = Spindeterminant.bitstring @@ Determinant.beta ki in\n",
    "            let j = Spindeterminant.bitstring @@ Determinant.beta kj in\n",
    "            Bitstring.to_list (Bitstring.logand i j)\n",
    "          in\n",
    "          match s with\n",
    "          | Spin.Alfa -> alfa, beta\n",
    "          | Spin.Beta -> beta, alfa\n",
    "        in\n",
    "\n",
    "        if s = s' then            \n",
    "            let s'' = Spin.other s' in\n",
    "            sum mos_cabs (fun a -> \n",
    "            sum mos  (fun n ->\n",
    "                 h_two i n a k s  s  *.  f_two j a n l s  s \n",
    "              +. h_two i n a l s  s  *.  f_two j a k n s  s \n",
    "              -. h_two j n a k s  s  *.  f_two i a n l s  s   \n",
    "              -. h_two j n a l s  s  *.  f_two i a k n s  s\n",
    "              )\n",
    "            +. sum mos_virt (fun m ->\n",
    "                -. h_two i j a m s  s  *.  f_two m a k l s  s )\n",
    "            +. sum mos' (fun n ->\n",
    "                 h_two i n k a s  s'' *.  f_two j a l n s  s''\n",
    "              +. h_two j n l a s  s'' *.  f_two i a k n s  s''\n",
    "              -. h_two i n l a s  s'' *.  f_two j a k n s  s''\n",
    "              -. h_two j n k a s  s'' *.  f_two i a l n s  s''\n",
    "              )\n",
    "            )\n",
    "        else\n",
    "            sum mos_cabs (fun a ->\n",
    "            sum mos_virt' (fun m ->\n",
    "                  h_two i j a m s  s' *.  f_two a m k l s  s' ) +.\n",
    "            sum mos_virt  (fun m ->\n",
    "                  h_two i j m a s  s' *.  f_two m a k l s  s' ) +.\n",
    "            sum mos  (fun n ->\n",
    "                  h_two n i a k s  s  *.  f_two a j n l s  s'\n",
    "               +. h_two n j a l s  s' *.  f_two i a k n s  s \n",
    "               -. h_two n j k a s  s' *.  f_two i a n l s  s'\n",
    "               ) +.\n",
    "            sum mos' (fun n -> if n >= j then 0. else \n",
    "                  h_two i n k a s  s' *.  f_two j a l n s' s'\n",
    "               +. h_two n j a l s' s' *.  f_two i a k n s  s' ) +. \n",
    "            sum mos' (fun n -> if n <= j then 0. else \n",
    "               -. h_two i n k a s  s' *.  f_two j a n l s' s' \n",
    "               -. h_two j n a l s' s' *.  f_two i a k n s  s' ) +.\n",
    "            sum mos' (fun n -> \n",
    "               -. h_two i n a l s  s' *.  f_two a j k n s  s'\n",
    "               )\n",
    "            )\n",
    "    in\n",
    "\n",
    "    match phase with\n",
    "        | Phase.Pos -> result\n",
    "        | Phase.Neg -> -. result\n",
    "\n",
    "\n",
    "let _ = check 300 integral_value 2 1 2 2"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/markdown",
    "scrolled": false
   },
   "source": [
    "let ki = det_I.(99)\n",
    "let kj = det_I.(193)\n",
    "\n",
    "let alpha_to_string alpha = \n",
    "   \n",
    "      let exc0 = Array.init (aux_num+1) (fun _ -> [|\"-\";\"-\"|]) in\n",
    " \n",
    "      let i, j, k, l, s, s', phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      let spin = function\n",
    "        | Spin.Alfa -> 0\n",
    "        | _ -> 1\n",
    "      in\n",
    "        exc0.(i).(spin s ) <- \"i\" \n",
    "      ; exc0.(j).(spin s') <- \"j\" \n",
    "      ; exc0.(k).(spin s ) <- \"k\" \n",
    "      ; exc0.(l).(spin s') <- \"l\" \n",
    "      ;\n",
    "      let s0, s0' = s, s' in\n",
    "    \n",
    "      let i, j, k, l, s, s', p1 =\n",
    "          match Excitation.of_det ki alpha with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      \n",
    "      if exc0.(i).(spin s ) =  \"-\" then exc0.(i).(spin s ) <- \"n\";\n",
    "      if exc0.(j).(spin s') =  \"-\" then exc0.(j).(spin s') <- \"n\";\n",
    "      if exc0.(k).(spin s ) =  \"-\" then exc0.(k).(spin s ) <- if k > mo_num then \"a\" else \"m\";\n",
    "      if exc0.(l).(spin s') =  \"-\" then exc0.(l).(spin s') <- if l > mo_num then \"a\" else \"m\";\n",
    "      \n",
    "      let string_h = \n",
    "          Printf.sprintf \"h_two %s %s %s %s %s %s *.  \" \n",
    "            exc0.(i).(spin s )\n",
    "            exc0.(j).(spin s')\n",
    "            exc0.(k).(spin s )\n",
    "            exc0.(l).(spin s')\n",
    "            (if s  = s0 then \"s \" else \"s'\")\n",
    "            (if s' = s0 then \"s \" else \"s'\")\n",
    "      in \n",
    "      \n",
    "      let i, j, k, l, s, s', p2 =\n",
    "          match Excitation.of_det alpha kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      \n",
    "      let string_f = \n",
    "          Printf.sprintf \"f_two %s %s %s %s %s %s\" \n",
    "            exc0.(i).(spin s )\n",
    "            exc0.(j).(spin s')\n",
    "            exc0.(k).(spin s )\n",
    "            exc0.(l).(spin s')\n",
    "            (if s  = s0 then \"s \" else \"s'\")\n",
    "            (if s' = s0 then \"s \" else \"s'\")\n",
    "      in\n",
    "(*\n",
    "        Format.printf \"|I> -> |a> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det ki alpha) string_h ;\n",
    "        Format.printf \"|a> -> |J> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det alpha kj) string_f ;\n",
    "        *)\n",
    "\n",
    "    \n",
    "        let sign = \n",
    "            if Phase.add p1 p2 = phase then \"+.\" else \"-.\"\n",
    "        in\n",
    "        sign ^ string_h ^ string_f\n",
    "    \n",
    "\n",
    "let alpha_debug alpha = \n",
    "   \n",
    "      let i, j, k, l, s, s', p1 =\n",
    "          match Excitation.of_det ki alpha with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      Printf.printf \"%d %d %d %d \" i j k l;\n",
    "      \n",
    "      let i, j, k, l, s, s', p2 =\n",
    "          match Excitation.of_det alpha kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      (*\n",
    "        Format.printf \"|I> -> |a> : %a | \\n@.\" Excitation.pp (Excitation.of_det ki alpha) ;\n",
    "        Format.printf \"|a> -> |J> : %a | \\n@.\" Excitation.pp (Excitation.of_det alpha kj) ;\n",
    "        *)\n",
    "      Printf.printf \"%d %d %d %d \\n%!\" i j k l\n",
    "\n",
    "      \n",
    "let strings = \n",
    "    Format.printf \"|I> -> |J> : %a |\\n@.\" Excitation.pp (Excitation.of_det ki kj) ;\n",
    "    generate_alphas ki kj 2 1 2 2\n",
    "    |> Array.of_list\n",
    "    |> Array.mapi (fun kk alpha -> alpha_to_string alpha)\n",
    "    |> Array.to_list\n",
    "    |> List.sort_uniq compare\n",
    "    |> Array.of_list\n",
    "\n",
    "let _ = Array.iteri (fun i x -> Printf.printf \"%d %s \\n%!\" i x) strings\n",
    "\n",
    "let _ =\n",
    "  let v =\n",
    "    let alphas =\n",
    "      generate_alphas ki kj 2 1 2 2\n",
    "      (*\n",
    "      |> List.filter (fun alpha ->\n",
    "        let x =  alpha_to_string alpha in\n",
    "           x = strings.(6)\n",
    "        )\n",
    "        *)\n",
    "    in\n",
    "    (*\n",
    "    List.iter alpha_debug alphas ;\n",
    "    Printf.printf \"\\n%!\";\n",
    "    *)\n",
    "    compute_HaaF ki alphas kj \n",
    "  in\n",
    "  let x = (integral_value ki kj) in\n",
    "  Printf.printf \"%20.8e %20.8e %20.8e\\n%!\" x v (v-. x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12. 2 2 22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let _ = 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",
    "           sum mos_cabs (fun b -> \n",
    "           sum mos_cabs (fun a -> \n",
    "               h_two h h' a b s s'  *. f_two a b p p' s s'\n",
    "           )) \n",
    "       else (* Alpha-Alpha / Beta-Beta *)\n",
    "           sum mos_cabs (fun b -> \n",
    "           sum mos_cabs (fun a ->  if b >= a then 0. else\n",
    "               h_two h h' a b s s'  *. f_two a b p p' s s'\n",
    "           )) \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": [
    "let _ = png_image \"3_1_22.png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$|J\\rangle = \\hat{T}_{ijm}^{kln}|I\\rangle$\n",
    "\n",
    "$\\sum_{a}\n",
    "    \\hat{T}_{ij}^{al} \\hat{T}_{am}^{kn}   \n",
    "  + \\hat{T}_{ij}^{ka} \\hat{T}_{am}^{ln}   \n",
    "  + \\hat{T}_{im}^{ka} \\hat{T}_{ja}^{ln}   \n",
    "  + \\hat{T}_{im}^{an} \\hat{T}_{aj}^{kl}   \n",
    "  + \\hat{T}_{jm}^{an} \\hat{T}_{ia}^{kl}  \n",
    "  + \\hat{T}_{jm}^{la} \\hat{T}_{ia}^{kn} \n",
    "$\n",
    "\n",
    "$$\n",
    "\\sum_{a} \\langle i j || a l \\rangle  \\left[ a m || k n \\right] +\n",
    "         \\langle i j || k a \\rangle  \\left[ a m || l n \\right] +\n",
    "         \\langle i m || k a \\rangle  \\left[ j a || l n \\right] +\n",
    "         \\langle i m || a n \\rangle  \\left[ a j || k l \\right] +\n",
    "         \\langle j m || a n \\rangle  \\left[ i a || k l \\right] +\n",
    "         \\langle j m || l a \\rangle  \\left[ i a || k n \\right] \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "let integral_value ki kj = \n",
    "    let i, j, m, k, l, n, 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",
    "    let result = \n",
    "       let open Spin in \n",
    "       match s1, s2, s3 with\n",
    "       | Alfa, Alfa, Alfa\n",
    "       | Beta, Beta, Beta -> \n",
    "           sum mos_cabs (fun a -> \n",
    "                h_two i j a k s1 s2 *.  f_two m a l n s3 s3 \n",
    "             +. h_two i j a n s1 s2 *.  f_two m a k l s3 s2 \n",
    "             +. h_two i m a l s1 s3 *.  f_two j a k n s2 s3 \n",
    "             +. h_two j m a k s2 s3 *.  f_two i a l n s1 s3 \n",
    "             +. h_two j m a n s2 s3 *.  f_two i a k l s1 s2 \n",
    "             -. h_two i j a l s1 s2 *.  f_two m a k n s3 s3 \n",
    "             -. h_two i m a k s1 s3 *.  f_two j a l n s2 s3 \n",
    "             -. h_two i m a n s1 s3 *.  f_two j a k l s2 s2 \n",
    "             -. h_two j m a l s2 s3 *.  f_two i a k n s1 s3 )\n",
    "       | Alfa, Alfa, Beta\n",
    "       | Beta, Beta, Alfa -> \n",
    "           sum mos_cabs (fun a -> \n",
    "                h_two i j a l s1 s2 *.  f_two a m k n s1 s3\n",
    "             +. h_two i m k a s1 s3 *.  f_two j a l n s2 s3\n",
    "             +. h_two j m a n s2 s3 *.  f_two i a k l s1 s2 \n",
    "             +. h_two j m l a s2 s3 *.  f_two i a k n s1 s3\n",
    "             -. h_two i j a k s1 s2 *.  f_two a m l n s1 s3\n",
    "             -. h_two i m a n s1 s3 *.  f_two j a k l s2 s2 \n",
    "             -. h_two i m l a s1 s3 *.  f_two j a k n s2 s3\n",
    "             -. h_two j m k a s2 s3 *.  f_two i a l n s1 s3\n",
    "           ) \n",
    "       | Alfa, Beta, Beta\n",
    "       | Beta, Alfa, Alfa -> \n",
    "           sum mos_cabs (fun a -> \n",
    "                h_two i j a l s1 s2 *.  f_two a m k n s1 s3 \n",
    "             +. h_two i m a n s1 s3 *.  f_two a j k l s1 s2 \n",
    "             +. h_two i m k a s1 s3 *.  f_two j a l n s2 s3 \n",
    "             +. h_two j m a n s2 s3 *.  f_two i a k l s1 s2 \n",
    "             -. h_two i j a n s1 s2 *.  f_two a m k l s1 s2 \n",
    "             -. h_two i j k a s1 s2 *.  f_two m a l n s2 s3 \n",
    "             -. h_two i m a l s1 s3 *.  f_two a j k n s1 s3 \n",
    "             -. h_two j m a l s2 s3 *.  f_two i a k n s1 s3 \n",
    "           ) \n",
    "       | Beta, Alfa, Beta\n",
    "       | Alfa, Beta, Alfa ->  assert false (*TODO *)\n",
    "    in\n",
    "    match phase with\n",
    "    | Phase.Pos -> result\n",
    "    | Phase.Neg -> -. result\n",
    "     \n",
    "\n",
    "let _ = check 200 integral_value 3 1 2 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "let ki = det_I.(129)\n",
    "let kj = det_I.(349)\n",
    "\n",
    "let alpha_to_string alpha = \n",
    "   \n",
    "      let exc0 = Array.init (aux_num+1) (fun _ -> [|\"-\";\"-\"|]) in\n",
    " \n",
    "      let i, j, m, k, l, n,  s, s', s'',  phase =\n",
    "          match Excitation.of_det ki kj with\n",
    "          | Excitation.(Triple (phase,\n",
    "              { hole ; particle ; spin },\n",
    "              {hole=hole' ; particle=particle' ; spin=spin' }, \n",
    "              {hole=hole''; particle=particle''; spin=spin''}  )) ->\n",
    "            hole, hole', hole'', particle, particle', particle'', spin, spin', spin'', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      let spin = function\n",
    "        | Spin.Alfa -> 0\n",
    "        | _ -> 1\n",
    "      in\n",
    "        exc0.(i).(spin s  ) <- \"i\" \n",
    "      ; exc0.(j).(spin s' ) <- \"j\" \n",
    "      ; exc0.(k).(spin s  ) <- \"k\" \n",
    "      ; exc0.(l).(spin s' ) <- \"l\" \n",
    "      ; exc0.(m).(spin s'') <- \"m\" \n",
    "      ; exc0.(n).(spin s'') <- \"n\" \n",
    "      ;\n",
    "      let s0, s0', s0'' = s, s', s'' in\n",
    "    \n",
    "      let i, j, k, l, s, s', p1 =\n",
    "          match Excitation.of_det ki alpha with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      \n",
    "      if exc0.(i).(spin s ) =  \"-\" then exc0.(i).(spin s ) <- \"p\";\n",
    "      if exc0.(j).(spin s') =  \"-\" then exc0.(j).(spin s') <- \"p\";\n",
    "      if exc0.(k).(spin s ) =  \"-\" then exc0.(k).(spin s ) <- if k > mo_num then \"a\" else \"q\";\n",
    "      if exc0.(l).(spin s') =  \"-\" then exc0.(l).(spin s') <- if l > mo_num then \"a\" else \"q\";\n",
    "      \n",
    "      let string_h = \n",
    "          Printf.sprintf \"h_two %s %s %s %s %s %s *.  \" \n",
    "            exc0.(i).(spin s )\n",
    "            exc0.(j).(spin s')\n",
    "            exc0.(k).(spin s )\n",
    "            exc0.(l).(spin s')\n",
    "            (if s  = s0  then \"s \" else if s = s0'  then \"s'\"  else \"s''\")\n",
    "            (if s' = s0' then \"s'\" else if s = s0'' then \"s''\" else \"s\"  )\n",
    "      in \n",
    "      \n",
    "      let i, j, k, l, s, s', p2 =\n",
    "          match Excitation.of_det alpha kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      \n",
    "      let string_f = \n",
    "          Printf.sprintf \"f_two %s %s %s %s %s %s\" \n",
    "            exc0.(i).(spin s )\n",
    "            exc0.(j).(spin s')\n",
    "            exc0.(k).(spin s )\n",
    "            exc0.(l).(spin s')\n",
    "            (if s  = s0  then \"s \" else if s = s0'  then \"s'\"  else \"s''\")\n",
    "            (if s' = s0' then \"s'\" else if s = s0'' then \"s''\" else \"s\"  )\n",
    "      in\n",
    "        (*\n",
    "        Format.printf \"|I> -> |a> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det ki alpha) string_h ;\n",
    "        Format.printf \"|a> -> |J> : %a | %s\\n@.\" Excitation.pp (Excitation.of_det alpha kj) string_f ;\n",
    "        *)\n",
    "\n",
    "    \n",
    "        let sign = \n",
    "            if Phase.add p1 p2 = phase then \"+.\" else \"-.\"\n",
    "        in\n",
    "        sign ^ string_h ^ string_f\n",
    "    \n",
    "\n",
    "let alpha_debug alpha = \n",
    "   \n",
    "      let i, j, k, l, s, s', p1 =\n",
    "          match Excitation.of_det ki alpha with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      Printf.printf \"%d %d %d %d \" i j k l;\n",
    "      \n",
    "      let i, j, k, l, s, s', p2 =\n",
    "          match Excitation.of_det alpha kj with\n",
    "          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->\n",
    "            hole, hole', particle, particle', spin, spin', phase\n",
    "          | _ -> assert false\n",
    "      in\n",
    "      (*\n",
    "        Format.printf \"|I> -> |a> : %a | \\n@.\" Excitation.pp (Excitation.of_det ki alpha) ;\n",
    "        Format.printf \"|a> -> |J> : %a | \\n@.\" Excitation.pp (Excitation.of_det alpha kj) ;\n",
    "        *)\n",
    "      Printf.printf \"%d %d %d %d \\n%!\" i j k l\n",
    "\n",
    "      \n",
    "let strings = \n",
    "    Format.printf \"|I> -> |J> : %a |\\n@.\" Excitation.pp (Excitation.of_det ki kj) ;\n",
    "    generate_alphas ki kj 3 1 2 2\n",
    "    |> Array.of_list\n",
    "    |> Array.mapi (fun kk alpha -> alpha_to_string alpha)\n",
    "    |> Array.to_list\n",
    "    |> List.sort_uniq compare\n",
    "    |> Array.of_list\n",
    "\n",
    "let _ = Array.iteri (fun i x -> Printf.printf \"%d %s \\n%!\" i x) strings\n",
    "\n",
    "let _ =\n",
    "  let v =\n",
    "    let alphas =\n",
    "      generate_alphas ki kj 3 1 2 2\n",
    "      (*\n",
    "      |> List.filter (fun alpha ->\n",
    "        let x =  alpha_to_string alpha in\n",
    "           x = strings.(6)\n",
    "        )\n",
    "        *)\n",
    "    in\n",
    "    (*\n",
    "    List.iter alpha_debug alphas ;\n",
    "    Printf.printf \"\\n%!\";\n",
    "    *)\n",
    "    compute_HaaF ki alphas kj \n",
    "  in\n",
    "  let x = (integral_value ki kj) in\n",
    "  Printf.printf \"%20.8e %20.8e %20.8e\\n%!\" x v (v-. x)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "OCaml 4.07.1",
   "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
}