diff --git a/Notebooks/F12_matrix.ipynb b/Notebooks/F12_matrix.ipynb
index 66ce4fe..41a40ee 100644
--- a/Notebooks/F12_matrix.ipynb
+++ b/Notebooks/F12_matrix.ipynb
@@ -1523,7 +1523,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# 8. 1 2 22a"
+    "# 8. 1 2 22"
    ]
   },
   {
@@ -1599,7 +1599,7 @@
     "    | Phase.Neg -> -. result\n",
     "     \n",
     "\n",
-    "let _ = check 50 integral_value 1 2 2 2"
+    "let _ = check 0 integral_value 1 2 2 2"
    ]
   },
   {
@@ -1688,13 +1688,345 @@
     "let _ = check 0 integral_value 2 1 1 2"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 10. 2 1 21"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "png_image \"2_1_21.png\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$|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}_{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 i j k a s s' *.\n",
+    "             ( f_one a l s' +.\n",
+    "               sum mos  (fun n -> f_two n a n l s' s') +.\n",
+    "               sum mos' (fun n -> f_two n a n l s'' s') \n",
+    "             ) \n",
+    "             ) \n",
+    "        in\n",
+    "\n",
+    "        match phase with\n",
+    "        | Phase.Pos -> result\n",
+    "        | Phase.Neg -> -. result\n",
+    "\n",
+    "\n",
+    "let _ = check 0 integral_value 2 1 2 1"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 11. *2 1 22*"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "png_image \"2_1_22.png\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$|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}_{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",
+    "$\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, mos_virt =\n",
+    "          match s with\n",
+    "          | Spin.Alfa -> mos_a ki, mos_virt_a\n",
+    "          | Spin.Beta -> mos_b ki, mos_virt_b\n",
+    "        in\n",
+    "        \n",
+    "        let mos', mos_virt' =\n",
+    "          match s' with\n",
+    "          | Spin.Alfa -> mos_a ki, mos_virt_a\n",
+    "          | Spin.Beta -> mos_b ki, mos_virt_b\n",
+    "        in\n",
+    "        \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 i n a l s s' *. f_two a j k n s s' ) -.\n",
+    "            sum mos (fun n -> \n",
+    "                 h_two n j k a s s' *. f_two i a n l s s' )\n",
+    "        )\n",
+    "    in\n",
+    "\n",
+    "    match phase with\n",
+    "        | Phase.Pos -> result\n",
+    "        | Phase.Neg -> -. result\n",
+    "\n",
+    "\n",
+    "let _ = check 0 integral_value 2 1 2 2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 12. 2 2 22"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "png_image \"2_2_22.png\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
+    "\n",
+    "$$\n",
+    "\\sum_{b}\\sum_{a<b}\n",
+    "    \\langle i j || a b \\rangle \n",
+    "    \\left[ ab || kj \\right]\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "let integral_value ki kj = \n",
+    "    let h, h', p, p', s, s', phase =\n",
+    "          match Excitation.of_det ki kj with\n",
+    "          | Excitation.(Double (phase,\n",
+    "              { hole=h ; particle=p ; spin=s },\n",
+    "              { hole=h'; particle=p'; spin=s'}) )  -> h, h', p, p', s, s', phase\n",
+    "          | _ -> assert false\n",
+    "    in\n",
+    "\n",
+    "    let result = \n",
+    "       if s <> s' then (* Alpha-Beta *)\n",
+    "           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": [
+    "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",
+    "       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 j k a s1 s2  *. f_two a m l n s2 s3 +.\n",
+    "           h_two i m k a s1 s3  *. f_two j a l n s2 s3 +.\n",
+    "           h_two i m a n s1 s3  *. f_two a j k l s1 s2 +.\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",
+    "           ) \n",
+    "    in\n",
+    "    match phase with\n",
+    "    | Phase.Pos -> result\n",
+    "    | Phase.Neg -> -. result\n",
+    "     \n",
+    "\n",
+    "let _ = check 0 integral_value 3 1 2 2"
+   ]
+  },
+  {
+   "cell_type": "raw",
    "metadata": {
     "scrolled": true
    },
-   "outputs": [],
    "source": [
     "let ki = det_I.(0)\n",
     "let kj = det_I.(11)\n",
@@ -1787,172 +2119,6 @@
     "\n"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# 10. *2 1 21*"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "png_image \"2_1_21.png\""
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# 11. *2 1 22*"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "png_image \"2_1_22.png\""
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# 12. 2 2 22"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "png_image \"2_2_22.png\""
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "$|J\\rangle = \\hat{T}_{ij}^{kl}|I\\rangle$\n",
-    "\n",
-    "$$\n",
-    "\\sum_{b}\\sum_{a<b}\n",
-    "    \\langle i j || a b \\rangle \n",
-    "    \\left[ ab || kj \\right]\n",
-    "$$"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "let integral_value ki kj = \n",
-    "    let h, h', p, p', s, s', phase =\n",
-    "          match Excitation.of_det ki kj with\n",
-    "          | Excitation.(Double (phase,\n",
-    "              { hole=h ; particle=p ; spin=s },\n",
-    "              { hole=h'; particle=p'; spin=s'}) )  -> h, h', p, p', s, s', phase\n",
-    "          | _ -> assert false\n",
-    "    in\n",
-    "\n",
-    "    let result = \n",
-    "       if s <> s' then (* Alpha-Beta *)\n",
-    "           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": [
-    "png_image \"3_1_22.png\""
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "$|J\\rangle = \\hat{T}_{ijm}^{kln}|I\\rangle$\n",
-    "\n",
-    "$$\n",
-    "\\sum_{a} \\langle i j || k a \\rangle  \\left[ m a || n l \\right] \n",
-    "$$"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "let integral_value ki kj = \n",
-    "    let h1, h2, h3, p1, p2, p3, s1, s2, s3, phase =\n",
-    "          match Excitation.of_det ki kj with\n",
-    "          | Excitation.(Triple (phase,\n",
-    "              { hole=h1 ; particle=p1 ; spin=s1 },\n",
-    "              { hole=h2 ; particle=p2 ; spin=s2 },\n",
-    "              { hole=h3 ; particle=p3 ; spin=s3 }) )  -> h1, h2, h3, p1, p2, p3, s1, s2, s3, phase\n",
-    "          | _ -> assert false\n",
-    "    in\n",
-    "\n",
-    "(*\n",
-    "123\n",
-    "213\n",
-    "231\n",
-    "321\n",
-    "312\n",
-    "132\n",
-    "*)\n",
-    "\n",
-    "    let result = \n",
-    "       List.fold_left (fun accu a -> accu +.\n",
-    "           h_two h1 h2 p1 a s1 s2  *. f_two h3 a p3 p2 s3 s2 -.\n",
-    "           h_two h2 h1 p2 a s2 s1  *. f_two h3 a p3 p1 s3 s1 +.\n",
-    "           h_two h2 h3 p2 a s2 s3  *. f_two h1 a p1 p3 s1 s3 -.\n",
-    "           h_two h3 h2 p3 a s3 s2  *. f_two h1 a p1 p2 s1 s2 +.\n",
-    "           h_two h3 h1 p3 a s3 s1  *. f_two h2 a p2 p1 s2 s1 -.\n",
-    "           h_two h1 h3 p1 a s1 s3  *. f_two h2 a p2 p3 s2 s3 \n",
-    "           ) 0. mos_cabs\n",
-    "    in\n",
-    "    match phase with\n",
-    "    | Phase.Pos -> result\n",
-    "    | Phase.Neg -> -. result\n",
-    "     \n",
-    "\n",
-    "let _ = check 100 integral_value 3 1 2 2"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,