diff --git a/Work.ipynb b/Work.ipynb
index 98917f4..922d530 100644
--- a/Work.ipynb
+++ b/Work.ipynb
@@ -78,7 +78,7 @@
     "(* --------- *)\n",
     "\n",
     "(*Mettre le bon chemin ici *)\n",
-    "#cd \"/home/scemama/QCaml\";;\n",
+    "#cd \"/home/ydamour/QCaml\";;\n",
     "\n",
     "#use \"topfind\";;\n",
     "#require \"jupyter.notebook\";;\n",
@@ -615,7 +615,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "let nocc = 2  (* On a 10 electrons, donc 5 orbitales occupees. *)\n",
+    "let nocc = 1  (* On a 10 electrons, donc 5 orbitales occupees. *)\n",
     "\n",
     "\n",
     "let c_of_h  m_H = \n",
@@ -1029,6 +1029,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "(*\n",
     "(* Test Boys *)\n",
     "\n",
     "let n_ao = Mat.dim1 m_C ;;\n",
@@ -1070,7 +1071,8 @@
     "        asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )));;\n",
     "\n",
     "\n",
-    "m_alpha_boys multipoles;;\n"
+    "m_alpha_boys multipoles;;\n",
+    "*)"
    ]
   },
   {
@@ -1244,6 +1246,24 @@
     "\n",
     "$= \\left(\\sum_e c_{ei} \\sum_f c_{fi}  - \\sum_e c_{ej} \\sum_f c_{fj} \\right) \\left(\\sum_a c_{ai} \\sum_b c_{bi}  - \\sum_a c_{aj} \\sum_b c_{bj} \\right) [ \\chi_a \\chi_b| \\chi_e \\chi_f ] $\n",
     "\n",
+    "(\n",
+    "Ainsi comme :\n",
+    "\n",
+    "$\\phi_i ^2 = \\left( \\sum_a c_{ai} \\chi_a \\right)^2$\n",
+    "$= \\sum_a c_{ai} \\sum_b c_{bi} \\chi_a \\chi_b$\n",
+    "\n",
+    "On aura :\n",
+    "\n",
+    "$D=\\sum_n [\\phi_n^2|\\phi_n^2]$\n",
+    "\n",
+    "$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} [\\chi_a \\chi_b| \\phi_n^2])$\n",
+    "\n",
+    "$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} \\sum_e c_{en} \\sum_f c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$\n",
+    "\n",
+    "$= \\sum_n (\\sum_a \\sum_b \\sum_e \\sum_f c_{an} c_{bn} c_{en} c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$\n",
+    "\n",
+    ")\n",
+    "\n",
     "Mais aussi :\n",
     "\n",
     "$B_{ij} = [\\phi_i ^2 - \\phi_j ^2 | \\phi_i \\phi_j ] = [\\left( \\sum_a c_{ai} \\chi_i \\right)^2 - \\left( \\sum_a c_{aj} \\chi_j \\right)^2| \\phi_i \\phi_j ] $\n",
@@ -1285,86 +1305,6 @@
     "## Calculs : Localisation de Edminston"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "(* Localisation de Edminstion *)\n",
-    "let ee_ints = AOBasis.ee_ints ao_basis;;\n",
-    "let m_C = MOBasis.mo_coef mo_basis;;\n",
-    "let n_ao = Mat.dim1 m_C ;;\n",
-    "let n_mo = n_ao - 1;;\n",
-    "let sum a = \n",
-    "  Array.fold_left (fun accu x -> accu +. x) 0. a\n",
-    ";;\n",
-    "\n",
-    "(*Calcul des intégrales*)  \n",
-    "let integral_general g i j =\n",
-    "Array.map (fun a ->\n",
-    "   let v = \n",
-    "    Array.map (fun b ->\n",
-    "        let u = \n",
-    "          Array.map (fun e ->\n",
-    "            let t = Array.map (fun f ->\n",
-    "            (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
-    "           ) (Util.array_range 1 n_ao)\n",
-    "           in sum t\n",
-    "        ) (Util.array_range 1 n_ao)\n",
-    "        in sum u\n",
-    "    ) (Util.array_range 1 n_ao)\n",
-    "   in  sum v\n",
-    ") (Util.array_range 1 n_ao)\n",
-    "|> sum"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "(* Calcul de toutes les intégrales B_12 -> Matrice *)\n",
-    "let m_b12 = Mat.init_cols n_ao n_ao (fun i j -> \n",
-    "integral_general (fun a b e f i j ->\n",
-    "                ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )  *. m_C.{e,i} *. m_C.{f,j}\n",
-    "    ) i j\n",
-    "        )\n",
-    "\n",
-    "(* Calcul de toutes les intégrales A_12 -> Matrice *)\n",
-    "let m_a12 = Mat.init_cols n_ao n_ao (fun i j ->\n",
-    "integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i}  *. m_C.{e,i} *. m_C.{f,j} \n",
-    "       -. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} ) \n",
-    "       *. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )\n",
-    "    ) i j\n",
-    "        )\n",
-    "\n",
-    "(* Calcul de tous les alpha -> Matrice *)\n",
-    "let m_alpha = Mat.init_cols n_ao n_ao ( fun i j ->\n",
-    "    asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n",
-    "    )\n",
-    "        )"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "(* Extraction du plus grand angle alpha*)\n",
-    "let ij_max_alpha = Mat.as_vec m_alpha\n",
-    "|> iamax;;\n",
-    "\n",
-    "(* Indices du plus grand angle alpha *)\n",
-    "let indice_i = ij_max_alpha mod n_ao;;\n",
-    "let indice_j = ij_max_alpha / n_ao + 1;;\n",
-    "\n",
-    "(* Valeur du plus grand angle alpha*)\n",
-    "let alpha = m_alpha.{indice_i,indice_j};;"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -1468,7 +1408,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Rotation des orbitales"
+    "### ER-V1"
    ]
   },
   {
@@ -1477,6 +1417,100 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "\n",
+    "(* Localisation de Edminstion *)\n",
+    "let ee_ints = AOBasis.ee_ints ao_basis;;\n",
+    "let m_C = MOBasis.mo_coef mo_basis;;\n",
+    "let n_ao = Mat.dim1 m_C ;;\n",
+    "let n_mo = Mat.dim2 m_C ;;\n",
+    "let sum a = \n",
+    "  Array.fold_left (fun accu x -> accu +. x) 0. a\n",
+    ";;\n",
+    "\n",
+    "(*Calcul des intégrales*)  \n",
+    "let integral_general g i j =\n",
+    "Array.map (fun a ->\n",
+    "   let v = \n",
+    "    Array.map (fun b ->\n",
+    "        let u = \n",
+    "          Array.map (fun e ->\n",
+    "            let t = Array.map (fun f ->\n",
+    "            (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
+    "           ) (Util.array_range 1 n_ao)\n",
+    "           in sum t\n",
+    "        ) (Util.array_range 1 n_ao)\n",
+    "        in sum u\n",
+    "    ) (Util.array_range 1 n_ao)\n",
+    "   in  sum v\n",
+    ") (Util.array_range 1 n_ao)\n",
+    "|> sum\n",
+    ";;\n",
+    "\n",
+    "(* Calcul de toutes les intégrales B_12 -> Matrice *)\n",
+    "let m_b12 = Mat.init_cols n_ao n_ao (fun i j -> \n",
+    "integral_general (fun a b e f i j ->\n",
+    "                ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )  *. m_C.{e,i} *. m_C.{f,j}\n",
+    "    ) i j\n",
+    "        )\n",
+    "\n",
+    "(* Calcul de toutes les intégrales A_12 -> Matrice *)\n",
+    "let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->\n",
+    "integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j}  *. m_C.{e,i} *. m_C.{f,j} \n",
+    "       -. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} ) \n",
+    "       *. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )\n",
+    "    ) i j\n",
+    "        )\n",
+    "\n",
+    "(* Calcul de tous les alpha -> Matrice *)\n",
+    "let m_alpha = Mat.init_cols n_mo n_mo ( fun i j ->\n",
+    "    asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n",
+    "    )\n",
+    "        )\n",
+    "\n",
+    "let s_D = \n",
+    "    let v_D = \n",
+    "        let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
+    "        integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i}  *. m_C.{e,i} *. m_C.{f,i} \n",
+    "           ) i j\n",
+    "               )\n",
+    "    in Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
+    "in  Vec.sum v_D   \n",
+    "\n",
+    "(*\n",
+    "(* Calcul de D *)\n",
+    "let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
+    "integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i}  *. m_C.{e,i} *. m_C.{f,i} \n",
+    "       ) i j\n",
+    "        )\n",
+    "        \n",
+    "let v_D = Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
+    "\n",
+    "let s_D = Vec.sum v_D\n",
+    "*)\n",
+    "\n",
+    " (*       \n",
+    "(* Calcul de D *)\n",
+    "let v_D = Vec.init n_ao (fun i  ->\n",
+    "integral_general (fun a b e f i  -> m_C.{a,i} *. m_C.{b,i}  *. m_C.{e,i} *. m_C.{f,i} \n",
+    "       ) i \n",
+    "        )\n",
+    "\n",
+    "*)\n",
+    "(* Rotation des orbitales *)\n",
+    "(*\n",
+    "(* Extraction du plus grand angle alpha*)\n",
+    "let ij_max_alpha = Mat.as_vec m_alpha\n",
+    "|> iamax;;\n",
+    "\n",
+    "(* Indices du plus grand angle alpha *)\n",
+    "let indice_i = ij_max_alpha mod n_ao;;\n",
+    "let indice_j = ij_max_alpha / n_ao + 1;;\n",
+    "\n",
+    "(* Valeur du plus grand angle alpha*)\n",
+    "let alpha = m_alpha.{indice_i,indice_j};;\n",
+    "*)\n",
+    "\n",
+    "(*\n",
     "let m_Phi = m_C\n",
     "\n",
     "(* Matrice de rotation 2 par 2 *)\n",
@@ -1492,22 +1526,12 @@
     "\n",
     "(* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle on\n",
     "obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n",
-    "let m_Ksi_tilde = gemm m_Ksi m_R;;"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Réinjection"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    "let m_Ksi_tilde = gemm m_Ksi m_R;;\n",
+    "*)\n",
+    "\n",
+    "(* Réinjection*)\n",
+    "\n",
+    "(*\n",
     "(* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur \n",
     "les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice\n",
     "des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec \n",
@@ -1528,7 +1552,8 @@
     "let m_interm = Mat.sub m_Phi m_Psi;;\n",
     "\n",
     "(* Construction de la nouvelle matrice d'OMs phi~ *)\n",
-    "let m_Phi_tilde = Mat.add m_Psi_tilde m_interm;; "
+    "let m_Phi_tilde = Mat.add m_Psi_tilde m_interm;; \n",
+    "*)"
    ]
   },
   {
@@ -1538,13 +1563,6 @@
     "# TESTS"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -1558,6 +1576,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "\n",
     "(* Fonction pour itérer , cc critère de convergence*) \n",
     "(* Localisation de Edminstion *)\n",
     "let ee_ints = AOBasis.ee_ints ao_basis;;\n",
@@ -1592,16 +1611,23 @@
     "  loc  : float ;\n",
     "}\n",
     "\n",
-    "let prev_iter = {coef = m_C ;loc = 0. }\n",
+    "(*let prev_iter = {coef = m_C ;loc = 0. }*)\n",
     "\n",
     "(* Calcul de la nouvelle matrice Phi après rotations *)\n",
     "let rec new_Phi prev_iter n = \n",
     " Printf.printf \"%d\\n%!\" n;\n",
+    " (*\n",
+    " Util.debug_matrix \"prev_iter\" prev_iter.coef;\n",
+    " *)\n",
     "    if n == 0 then\n",
     "      prev_iter\n",
     "    else\n",
     "    let m_Phi = prev_iter.coef\n",
     "    in\n",
+    "    (*\n",
+    "    Util.debug_matrix \"m_coef\" prev_iter.coef;\n",
+    "    Util.debug_matrix \"m_Phi\" m_Phi;\n",
+    "    *)\n",
     "    (* Calcul de tous les alpha -> Matrice *)\n",
     "    let f_alpha m_Phi =\n",
     "\n",
@@ -1614,49 +1640,70 @@
     "        in\n",
     "        (* Calcul de toutes les intégrales A_12 -> Matrice *)\n",
     "        let m_a12  = Mat.init_cols n_mo n_mo (fun i j ->\n",
-    "            integral_general (fun a b e f i j -> m_Phi.{a,i} *. m_Phi.{b,i}  *. m_Phi.{e,i} *. m_Phi.{f,j} \n",
+    "            integral_general (fun a b e f i j -> m_Phi.{a,i} *. m_Phi.{b,j}  *. m_Phi.{e,i} *. m_Phi.{f,j} \n",
     "                -. 0.25 *. ( m_Phi.{e,i} *. m_Phi.{f,i} -. m_Phi.{e,j} *. m_Phi.{f,j} ) \n",
     "               *. ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} )\n",
     "                ) i j\n",
     "                    )\n",
     "        in\n",
-    "        Mat.init_cols n_ao n_mo ( fun i j ->\n",
+    "        Mat.init_cols n_mo n_mo ( fun i j ->\n",
     "           asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n",
     "                )\n",
     "                    )\n",
-    "\n",
-    "    in\n",
+    "        in \n",
+    "        (* Calcul de D *)\n",
+    "        let s_D = \n",
+    "            let v_D = \n",
+    "                let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
+    "                integral_general (fun a b e f i j -> m_Phi.{a,i} *. m_Phi.{b,i}  *. m_Phi.{e,i} *. m_Phi.{f,i} \n",
+    "                   ) i j\n",
+    "                       )\n",
+    "            in Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
+    "        in  Vec.sum v_D \n",
+    "                      \n",
+    "    in  \n",
     "    (* Matrice de rotation 2 par 2 *)\n",
-    "    let indice_i, indice_j, m_R = \n",
+    "    let indice_i, indice_j, m_R  = \n",
     "\n",
+    "        \n",
     "        let m_alpha = f_alpha m_Phi\n",
     "        in\n",
+    "        (*\n",
+    "        Util.debug_matrix \"m_alpha\" m_alpha;\n",
+    "        *)\n",
     "        (* Valeur du plus grand angle alpha*)\n",
-    "        let alpha = \n",
-    "\n",
-    "            (* Extraction du plus grand angle alpha*)\n",
+    "        (* Extraction du plus grand angle alpha*)\n",
     "            let ij_max_alpha = Mat.as_vec m_alpha\n",
     "            |> iamax\n",
     "            in\n",
     "\n",
     "            (* Indices du plus grand angle alpha *)\n",
-    "            let indice_i = ij_max_alpha mod n_ao\n",
+    "            let indice_i = ij_max_alpha mod n_mo\n",
     "            in\n",
-    "            let indice_j = ij_max_alpha / n_ao + 1\n",
+    "            let indice_j = ij_max_alpha / n_mo + 1\n",
     "\n",
     "            in\n",
-    "            m_alpha.{indice_i,indice_j}\n",
+    "            (*\n",
+    "            Printf.printf \"%i\\n%!\" indice_i;\n",
+    "            Printf.printf \"%i\\n%!\" indice_j;\n",
+    "            *)\n",
+    "        let alpha = \n",
+    "           m_alpha.{indice_i,indice_j}\n",
     "        in\n",
+    "        \n",
+    "        (*\n",
     "        Printf.printf \"%f\\n%!\" alpha;\n",
-    "    \n",
+    "        *)\n",
     "        indice_i, indice_j, \n",
     "        Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha\n",
     "                                                else if i>j then sin alpha \n",
     "                                                else -. sin alpha )\n",
     "\n",
     "    in\n",
+    "    (*\n",
     "    Util.debug_matrix \"m_R\" m_R;\n",
     "    Printf.printf \"%d %d\\n%!\" indice_i indice_j;\n",
+    "    *)\n",
     "    (* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle \n",
     "    on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n",
     "    let m_Ksi, m_Ksi_tilde = \n",
@@ -1672,9 +1719,10 @@
     "        m_Ksi, gemm m_Ksi m_R\n",
     "\n",
     "    in\n",
-    "    \n",
+    "    (*\n",
     "    Util.debug_matrix \"m_Ksi\" m_Ksi;\n",
     "    Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n",
+    "    *)\n",
     "    (* Construction de la nouvelle matrice d'OMs phi~ *)\n",
     "    let m_Phi_tilde m_Ksi m_Ksi_tilde = \n",
     "\n",
@@ -1700,20 +1748,31 @@
     "    in    \n",
     "    let coef = m_Phi_tilde m_Ksi m_Ksi_tilde\n",
     "    in\n",
+    "    \n",
     "    Util.debug_matrix \"coef\" coef;\n",
-    "    let loc =  alpha \n",
+    "    \n",
+    "    let loc =  s_D \n",
     "    in\n",
-    "    let result = \n",
+    "    \n",
+    "    Printf.printf \"%f\\n%!\" loc;\n",
+    "    Printf.printf \"%f\\n%!\" prev_iter.loc;\n",
+    "    \n",
     "    if abs_float (loc -. prev_iter.loc)  < 1.e-5 \n",
     "        then\n",
     "            { coef  ; loc }\n",
     "        else\n",
-    "            prev_iter\n",
-    "    in\n",
-    "    new_Phi result (n-1)\n",
-    "\n",
-    "\n",
+    "            new_Phi {coef; loc;} (n-1)\n",
     "\n",
+    "    \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "new_Phi {coef = m_C ;loc = 0. } 10;;\n",
     "\n"
    ]
   },
@@ -1723,7 +1782,199 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    " new_Phi prev_iter 3;;"
+    "            "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "(* Fonction pour itérer , cc critère de convergence*) \n",
+    "(* Localisation de Edminstion *)\n",
+    "let ee_ints = AOBasis.ee_ints ao_basis;;\n",
+    "let m_C = MOBasis.mo_coef mo_basis;;\n",
+    "let n_ao = Mat.dim1 m_C ;;\n",
+    "let n_mo = Mat.dim2 m_C ;;\n",
+    "let sum a = \n",
+    "  Array.fold_left (fun accu x -> accu +. x) 0. a\n",
+    ";;\n",
+    "\n",
+    "(*Calcul des intégrales*)  \n",
+    "let integral_general g i j =\n",
+    "Array.map (fun a ->\n",
+    "   let v = \n",
+    "    Array.map (fun b ->\n",
+    "        let u = \n",
+    "          Array.map (fun e ->\n",
+    "            let t = Array.map (fun f ->\n",
+    "            (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
+    "           ) (Util.array_range 1 n_ao)\n",
+    "           in sum t\n",
+    "        ) (Util.array_range 1 n_ao)\n",
+    "        in sum u\n",
+    "    ) (Util.array_range 1 n_ao)\n",
+    "   in  sum v\n",
+    ") (Util.array_range 1 n_ao)\n",
+    "|> sum   \n",
+    "\n",
+    "let m_Phi = None;;\n",
+    "type iteration = \n",
+    "{ coef : Mat.t ;\n",
+    "  loc  : float ;\n",
+    "}\n",
+    "\n",
+    "(*let prev_iter = {coef = m_C ;loc = 0. }*)\n",
+    "\n",
+    "(* Calcul de la nouvelle matrice Phi après rotations *)\n",
+    "let rec new_Phi prev_iter n = \n",
+    " Printf.printf \"%d\\n%!\" n;\n",
+    " (*\n",
+    " Util.debug_matrix \"prev_iter\" prev_iter.coef;\n",
+    " *)\n",
+    "    if n == 0 then\n",
+    "      prev_iter\n",
+    "    else\n",
+    "    let m_Phi = prev_iter.coef\n",
+    "    in\n",
+    "    (*\n",
+    "    Util.debug_matrix \"m_coef\" prev_iter.coef;\n",
+    "    Util.debug_matrix \"m_Phi\" m_Phi;\n",
+    "    *)\n",
+    "    (* Calcul de tous les alpha -> Matrice *)\n",
+    "    let f_alpha m_Phi =\n",
+    "\n",
+    "        (* Calcul de toutes les intégrales B_12 -> Matrice *)\n",
+    "        let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> \n",
+    "        integral_general (fun a b e f i j ->\n",
+    "            ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} )  *. m_Phi.{e,i} *. m_Phi.{f,j}\n",
+    "                ) i j\n",
+    "                    )\n",
+    "        in\n",
+    "        (* Calcul de toutes les intégrales A_12 -> Matrice *)\n",
+    "        let m_a12  = Mat.init_cols n_mo n_mo (fun i j ->\n",
+    "            integral_general (fun a b e f i j -> m_Phi.{a,i} *. m_Phi.{b,i}  *. m_Phi.{e,i} *. m_Phi.{f,j} \n",
+    "                -. 0.25 *. ( m_Phi.{e,i} *. m_Phi.{f,i} -. m_Phi.{e,j} *. m_Phi.{f,j} ) \n",
+    "               *. ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} )\n",
+    "                ) i j\n",
+    "                    )\n",
+    "        in\n",
+    "        Mat.init_cols n_mo n_mo ( fun i j ->\n",
+    "           asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n",
+    "                )\n",
+    "                    )\n",
+    "    in  \n",
+    "    (* Matrice de rotation 2 par 2 *)\n",
+    "    let indice_i, indice_j, m_R  = \n",
+    "        let n = 3\n",
+    "        in\n",
+    "        let rec new_alpha m_alpha n =\n",
+    "        \n",
+    "            let m_alpha = f_alpha m_Phi\n",
+    "            in\n",
+    "            (*\n",
+    "            Util.debug_matrix \"m_alpha\" m_alpha;\n",
+    "            *)\n",
+    "            (* Valeur du plus grand angle alpha*)\n",
+    "            (* Extraction du plus grand angle alpha*)\n",
+    "                let ij_max_alpha = Mat.as_vec m_alpha\n",
+    "                |> iamax\n",
+    "                in\n",
+    "\n",
+    "                (* Indices du plus grand angle alpha *)\n",
+    "                let indice_i = ij_max_alpha mod n_mo\n",
+    "                in\n",
+    "                let indice_j = ij_max_alpha / n_mo + 1\n",
+    "\n",
+    "                \n",
+    "                (*\n",
+    "                Printf.printf \"%i\\n%!\" indice_i;\n",
+    "                Printf.printf \"%i\\n%!\" indice_j;\n",
+    "                *)\n",
+    "            in\n",
+    "            if m_alpha.{indice_i,indice_j} < (3.14 /. 2.)\n",
+    "                then m_alpha\n",
+    "                else new_alpha (Mat.sub m_alpha (Mat.make n_ao n_ao (3.14 /. 2.)))\n",
+    "        \n",
+    "        in\n",
+    "        let alpha = m_alpha.{indice_i, indice_j}\n",
+    "         \n",
+    "        \n",
+    "        (*\n",
+    "        Printf.printf \"%f\\n%!\" alpha;\n",
+    "        *)\n",
+    "        indice_i, indice_j, \n",
+    "        Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha\n",
+    "                                                else if i>j then sin alpha \n",
+    "                                                else -. sin alpha )\n",
+    "\n",
+    "    in\n",
+    "    (*\n",
+    "    Util.debug_matrix \"m_R\" m_R;\n",
+    "    Printf.printf \"%d %d\\n%!\" indice_i indice_j;\n",
+    "    *)\n",
+    "    (* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle \n",
+    "    on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n",
+    "    let m_Ksi, m_Ksi_tilde = \n",
+    "\n",
+    "        (* Fonction d'extraction des 2 vecteurs propres i et j de la matrice des OMs pour les mettres dans la matrice Ksi \n",
+    "        (n par 2) pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n",
+    "        let p = indice_i\n",
+    "        in\n",
+    "        let q = indice_j\n",
+    "        in\n",
+    "        let m_Ksi = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_Phi.{i,p} else m_Phi.{i,q} )\n",
+    "        in\n",
+    "        m_Ksi, gemm m_Ksi m_R\n",
+    "\n",
+    "    in\n",
+    "    (*\n",
+    "    Util.debug_matrix \"m_Ksi\" m_Ksi;\n",
+    "    Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n",
+    "    *)\n",
+    "    (* Construction de la nouvelle matrice d'OMs phi~ *)\n",
+    "    let m_Phi_tilde m_Ksi m_Ksi_tilde = \n",
+    "\n",
+    "    \n",
+    "     (* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n",
+    "      let m_Psi_tilde = Mat.init_cols n_ao n_mo (fun i j -> if j=q then m_Ksi_tilde.{i,2}\n",
+    "                                                       else if j=p then m_Ksi_tilde.{i,1}\n",
+    "                                                        else 0.)\n",
+    "         in\n",
+    "        (* Matrice intérmédiaire où les orbitales i et j ont été supprimées et remplacées par des 0, par soustraction de la \n",
+    "        matrice Phi par la matrice *)\n",
+    "        let m_interm = \n",
+    "\n",
+    "            (* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *)                                            \n",
+    "            let m_Psi = Mat.init_cols n_ao n_mo (fun i j -> if j=q then m_Ksi.{i,2}\n",
+    "                                                        else if j=p then m_Ksi.{i,1}\n",
+    "                                                        else 0.)\n",
+    "\n",
+    "            in\n",
+    "            Mat.sub m_Phi m_Psi\n",
+    "        in\n",
+    "        Mat.add m_Psi_tilde m_interm\n",
+    "    in    \n",
+    "    let coef = m_Phi_tilde m_Ksi m_Ksi_tilde\n",
+    "    in\n",
+    "    \n",
+    "    Util.debug_matrix \"coef\" coef;\n",
+    "    \n",
+    "    let loc =  alpha \n",
+    "    in\n",
+    "    (*\n",
+    "    Printf.printf \"%f\\n%!\" loc;\n",
+    "    Printf.printf \"%f\\n%!\" prev_iter.loc;\n",
+    "    *)\n",
+    "    if abs_float (loc -. prev_iter.loc)  < 1.e-10 \n",
+    "        then\n",
+    "            { coef  ; loc }\n",
+    "        else\n",
+    "            new_Phi {coef; loc;} (n-1)\n",
+    "\n",
+    "               "
    ]
   },
   {