diff --git a/Basis/Multipole.ml b/Basis/Multipole.ml
index cbc9d44..630e211 100644
--- a/Basis/Multipole.ml
+++ b/Basis/Multipole.ml
@@ -24,6 +24,15 @@ let matrix_z2 t =  t.(5)
 let matrix_xy t =  t.(6)
 let matrix_yz t =  t.(7)
 let matrix_zx t =  t.(8)
+let matrix_x3 t =  t.(9)
+let matrix_y3 t =  t.(10)
+let matrix_z3 t =  t.(11)
+let matrix_x4 t =  t.(12)
+let matrix_y4 t =  t.(13)
+let matrix_z4 t =  t.(14)
+
+
+
 
 let cutoff = integrals_cutoff
 
@@ -37,7 +46,7 @@ let to_powers x =
 let contracted_class shell_a shell_b : float Zmap.t array =
 
   match Csp.make shell_a shell_b with
-  | None -> Array.init 9 (fun _ -> Zmap.create 0)
+  | None -> Array.init 15 (fun _ -> Zmap.create 0)
   | Some shell_p ->
     begin
 
@@ -45,7 +54,7 @@ let contracted_class shell_a shell_b : float Zmap.t array =
       let class_indices = Csp.zkey_array shell_p in
 
       let contracted_class =
-        Array.init 9 (fun _ -> Array.make (Array.length class_indices) 0.)
+        Array.init 15 (fun _ -> Array.make (Array.length class_indices) 0.)
       in
 
       let a_minus_b =
@@ -97,10 +106,26 @@ let contracted_class shell_a shell_b : float Zmap.t array =
                     Overlap_primitives.hvrr (Po.get xyz angMomA + 2, Po.get xyz angMomB)
                       expo_inv
                       (Co.get xyz a_minus_b, Co.get xyz center_pa)
+                  
+                   in
+                  (* 1D <i|(x-Xa)^3|j> *)
+                  let j k =
+                    let xyz = xyz_of_int k in
+                    Overlap_primitives.hvrr (Po.get xyz angMomA + 3, Po.get xyz angMomB)
+                      expo_inv
+                      (Co.get xyz a_minus_b, Co.get xyz center_pa)
+                   in
+                  (* 1D <i|(x-Xa)^4|j> *)
+                  let l k =
+                    let xyz = xyz_of_int k in
+                    Overlap_primitives.hvrr (Po.get xyz angMomA + 4, Po.get xyz angMomB)
+                      expo_inv
+                      (Co.get xyz a_minus_b, Co.get xyz center_pa)
+                   
                   in
                   let norm =  norm_coef_scales.(i) in
-                  let f0, f1, f2, g0, g1, g2, h0, h1, h2 =
-                    f 0, f 1, f 2, g 0, g 1, g 2, h 0, h 1, h 2
+                  let f0, f1, f2, g0, g1, g2, h0, h1, h2, j0, j1, j2 , l0, l1, l2 =
+                    f 0, f 1, f 2, g 0, g 1, g 2, h 0, h 1, h 2, j 0, j 1, j 2, l 0, l 1, l 2
                   in
                   let  x = g0 +. f0 *. xa in
                   let  y = g1 +. f1 *. ya in
@@ -108,6 +133,13 @@ let contracted_class shell_a shell_b : float Zmap.t array =
                   let x2 = h0 +. xa *. (2. *. x -. xa *. f0) in
                   let y2 = h1 +. ya *. (2. *. y -. ya *. f1) in
                   let z2 = h2 +. za *. (2. *. z -. za *. f2) in
+                  let x3 = j0 +. xa *. f0 *. (3. *. x2 -. 3. *. x  *. xa +. xa *. xa) in
+                  let y3 = j1 +. ya *. f1 *. (3. *. y2 -. 3. *. y  *. ya +. ya *. ya) in
+                  let z3 = j2 +. za *. f2 *. (3. *. z2 -. 3. *. z  *. za +. za *. za) in
+                  let x4 = l0 +. xa *. f0 *. ( 4. *. x3 -. 6. *. x2 *. xa +. 4. *. x *. xa *. xa -. xa *. xa *. xa) in
+                  let y4 = l1 +. ya *. f1 *. ( 4. *. y3 -. 6. *. y2 *. ya +. 4. *. y *. ya *. ya -. ya *. ya *. ya) in
+                  let z4 = l2 +. za *. f2 *. ( 4. *. z3 -. 6. *. z2 *. za +. 4. *. z *. za *. za -. za *. za *. za) in
+
                   let c = contracted_class in
                   let d = coef_prod *. norm in
                   c.(0).(i) <- c.(0).(i) +. d *.  x *. f1 *. f2;
@@ -119,6 +151,13 @@ let contracted_class shell_a shell_b : float Zmap.t array =
                   c.(6).(i) <- c.(6).(i) +. d *.  x *.  y *. f2;
                   c.(7).(i) <- c.(7).(i) +. d *. f0 *.  y *.  z;
                   c.(8).(i) <- c.(8).(i) +. d *.  x *. f1 *.  z;
+                  c.(9).(i) <- c.(9).(i) +. d *.  x3 *.  f1 *. f2;
+                  c.(10).(i) <- c.(10).(i) +. d *. f0 *. y3 *. f2;
+                  c.(11).(i) <- c.(11).(i) +. d *. f0 *. f1 *. z3;
+                  c.(12).(i) <- c.(12).(i) +. d *. x4 *. f1 *. f2;
+                  c.(13).(i) <- c.(13).(i) +. d *. f0 *. y4 *. f2;
+                  c.(14).(i) <- c.(14).(i) +. d *. f0 *. f1 *. z4;
+
                 ) class_indices
             end
         ) (Csp.coefs_and_shell_pairs shell_p);
@@ -147,7 +186,7 @@ let of_basis basis =
   and shell = Bs.contracted_shells basis
   in
 
-  let result = Array.init 9 (fun _ -> Mat.create n n) in
+  let result = Array.init 15 (fun _ -> Mat.create n n) in
   for j=0 to (Array.length shell) - 1 do
     for i=0 to j do
       (* Compute all the integrals of the class *)
@@ -155,7 +194,7 @@ let of_basis basis =
         contracted_class shell.(i) shell.(j)
       in
 
-      for k=0 to 8 do
+      for k=0 to 14 do
         Array.iteri (fun j_c powers_j ->
             let j_c = Cs.index shell.(j) + j_c + 1 in
             let xj = to_powers powers_j in
@@ -176,7 +215,7 @@ let of_basis basis =
       done;
     done;
   done;
-  for k=0 to 8 do
+  for k=0 to 14 do
     Mat.detri result.(k);
   done;
   result
diff --git a/Basis/Multipole.mli b/Basis/Multipole.mli
index 64d09de..74346bd 100644
--- a/Basis/Multipole.mli
+++ b/Basis/Multipole.mli
@@ -31,4 +31,22 @@ val matrix_y2 : t -> Mat.t
 val matrix_z2 : t -> Mat.t
 (** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *)
 
+val matrix_x3 : t -> Mat.t
+(** {% $$ \langle \chi_i | x^3 | \chi_j \rangle $$ %} *)
+
+val matrix_y3 : t -> Mat.t
+(** {% $$ \langle \chi_i | y^3 | \chi_j \rangle $$ %} *)
+
+val matrix_z3 : t -> Mat.t
+(** {% $$ \langle \chi_i | z^3 | \chi_j \rangle $$ %} *)
+
+val matrix_x4 : t -> Mat.t
+(** {% $$ \langle \chi_i | x^4 | \chi_j \rangle $$ %} *)
+
+val matrix_y4 : t -> Mat.t
+(** {% $$ \langle \chi_i | y^4 | \chi_j \rangle $$ %} *)
+
+val matrix_z4 : t -> Mat.t
+(** {% $$ \langle \chi_i | z^4 | \chi_j \rangle $$ %} *)
+
 val of_basis : Basis.t -> t
diff --git a/MOBasis/Localisation.ml b/MOBasis/Localisation.ml
new file mode 100644
index 0000000..e7d1148
--- /dev/null
+++ b/MOBasis/Localisation.ml
@@ -0,0 +1,1376 @@
+(* Algorithme de localisation des orbitales *)
+
+(* SOMMAIRE *)
+
+(* 0. Paramètres du calcul *)
+
+(* I. Définition des types *)
+
+(* II : Définition des fonctions pour la localisation *)
+
+    (* 1. Définitions de base nécessaire pour la suite *)
+    (* 2. Fonction de calcul de la matrice des angles selon la méthode de Edminston *)
+    (* 3. Fonction calcul de la matrice des angles selon la méthode de Boys façon Edminston  = NON FONCTIONNELLE  *)
+    (* 4. Fonction de calcul de la matrice des angles selon la méthode de Boys *)
+    (* 5. Pattern matching :  calcul de alpha et du critere D *)
+    (* 6. Norme d'une matrice *)
+    (* 7.  Détermination de alpha_max et ses indices i et j.
+           Si alpha max > pi/2 on soustrait pi/2 à la matrice des alphas de manière récursive
+           Si alpha max < pi/2 on ajoute pi/2 à la matrice des alphas de manière récursive *)
+    (* 8. Matrice de rotation 2 par 2 *)
+    (* 9. Matrice de rotation n par n avec un angle fixe  *)
+    (* 10. Fonction d'extraction des 2 vecteurs propres i et j de la matrice des OMs pour les mettres dans la matrice Ksi (n par 2)
+           pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)
+    (* 11. Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle
+           on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)
+    (* 12. Fonction pour la création de matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)
+    (* 13. Matrice intérmédiaire où les orbitales i et j ont été supprimées et remplacées par des 0, par soustraction de la matrice Phi
+           par la matrice *)
+
+(* III. Localisation des OMs *)
+
+    (* 1. Boucle de calcul de la nouvelle matrice des coefficient après n rotation d'orbitales  *)
+    (* 2. Fonction de localisation *)
+
+(* IV : Fonctions pour la séparation et le rassemblement des matrices *)
+
+    (* 1. Fonction de création d'une liste d'entier à partir d'un vecteur de float*)
+    (* 2. Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)
+    (* 3. Fonction de séparation d'une matrice en 2 sous matrice, la première matrice correspondant aux colonnes de la matrice dont le numéro est présent
+          dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)
+    (* 4. Liste des OMs occupées *)
+    (* 5. Liste des OMs de coeur *)
+    (* 6. Fonction de creation d'une liste des OMs pi virtuelles *)
+    (* 7. Fonction de rassemblement de 2 matrices *)
+
+(* V. Calcul *)
+
+    (* 1. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)
+    (* 1.1 Liste des OMs pi virtuelles *)
+    (* 2. Application d'une légère rotation d'un angle fixé aux OMs pour améliorer la convergence *)
+    (* 3. Légende *)
+    (* 4. Rotation des orbitales : localise localisation new_X / "methode" / epsilon = 1. si <1 diminue l'angle des rotations / nb itération / critère de convergence *)
+    (* 4.1 Fonction de localisation de toutes les orbitales virtuelles sans separation *)
+    (* 4.2 Fonction de localisation en separant les orbitales pi et sigma *)
+    (* 4.3 Fonction de localisation en separant les orbitales pi et sigma et les 10 dernières sigma *)
+    (* 4.4 Fonction séparant les orbitales pi, la premiere et la seconde moitié des orbitales sigma  *)
+    (* 5. Rassemblement des occupées et des virtuelles *)
+    (* 5.1 Rassemblement des orbitales occupées (coeur et non coeur) *)
+    (* 5.2 Pattern matching du type de separation des orbitales virtuelles et rassemblement des orbitales occupées et virtuelles *)
+
+(* VI. Analyse *)
+
+    (* 1. Fonction de calcul de l'extension spatiale des orbitales *)
+    (* 2. Fonction de tri par ordre croissant d'une liste *)
+    (* 3. Fonction de tri par ordre croissant de l'extension spatiale des orbitales *)
+    (* 4. Moyenne/variance/ecart-type/mediane de l'extension spatiale des orbitales *)
+    (* 5. Fonction d'affichage *)
+
+let methode = "Boys_er";; (* Method for th orbitals localization *)
+let iteration = 10000;;   (* Maximum number of iteration for each localization *)
+let cc = 10e-6;;          (* Convergence criterion *)
+let pre_angle = 0.001;; (* Rotation of the matrix with a small angle to improve the convergence *)
+let epsilon = 1.;;        (* Rotation angle = angle /. epsilon, default value : epsilon = 1. *)
+
+let x_end = 1;; (* Choice of the x_end for the pi/sigma/end separation, size of the second matrix with the sigma orbitals *)
+(* !!! WARNING : even if you don't use this localization method, 0 < x_end < total number of sigma orbitals.
+                 If you don't use this method, use the default value, x_end = 1   WARNING !!! *)
+
+(* Choice of the separation type of orbitals *)
+
+(* For valence ones *)
+let type_separation_occu = "pi/sigma";;      (* Separation of the pi and sigma orbitals *)
+(*let type_separation_occu = "1/2";;*)           (* Separation of the first and second half of the orbitals *)    
+(*let type_separation_occu = "pi/sigma/end";;*)  (* Separation pi,(n-x_end) sigma and x_end sigma *)
+(*let type_separation_occu = "pi/1/2";;*)        (* Separation pi, (n/2) first sigma and (n/2) last sigma *)
+(*let type_separation_occu = "full";;*)          (* Without separation *)
+
+(* For occupied ones *)
+let type_separation_vir = "pi/sigma";;         (* Separation of the pi and sigma orbitals *)
+(*let type_separation_vir = "1/2";;*)              (* Separation of the first and second half of the orbitals *)
+(*let type_separation_vir = "pi/sigma/end";;*)     (* Separation pi,(n-x_end) sigma and x_end sigma *)
+(*let type_separation_vir = "pi/1/2";;*)           (* Separation pi, (n/2) first sigma and (n/2) last sigma *)
+(*let type_separation_vir = "full";;*)             (* Without separation *)
+
+open Lacaml.D
+
+(* I. Types definitions *)
+
+type alphaij = {
+  alpha_max : float;
+  indice_ii : int;
+  indice_jj : int;
+}
+
+type alphad = {
+  m_alpha : Mat.t;
+  d : float;
+}
+
+let localize mo_basis =
+
+  let simulation = MOBasis.simulation mo_basis in
+  let nocc =
+      let elec = Simulation.electrons simulation
+          in
+          Electrons.n_alfa elec
+  in
+  let ao_basis = MOBasis.ao_basis mo_basis in
+  let n_core =
+    (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
+  in
+  let charge = Simulation.charge simulation
+               |>Charge.to_int  
+  in
+
+  (* II : Definition of the fonctions for the localization *)
+  
+  (* 1. Basic definitions *)
+  let m_C = MOBasis.mo_coef mo_basis in
+  let n_ao = Mat.dim1 m_C in
+  let pi = acos(-1.) in
+
+  (* 2. Function for Edminston-Ruedenberg localization*)
+  (* Coefficient matrix n_ao x n_mo -> n_mo x n_mo matrix with angle(rad) between orbitals, Edminston localization criterion *)
+  (* matrix -> matrix, float *)
+  let f_alpha_er m_C =
+      
+      let ee_ints = AOBasis.ee_ints ao_basis in
+      let n_mo = Mat.dim2 m_C in
+   
+      let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in
+      let m_a12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in
+      let v_d = Vec.init n_mo (fun i -> 0.) in
+   
+      (* Tableaux temporaires *)
+      let m_pqr =
+          Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)
+      in
+      let m_qr_i = Mat.create (n_ao*n_ao) n_mo in
+      let m_ri_j = Mat.create (n_ao*n_mo) n_mo in
+      let m_ij_k = Mat.create (n_mo*n_mo) n_mo in
+   
+      Array.iter (fun s ->
+         (* Grosse boucle externe sur s *)
+         Array.iter (fun r ->
+           Array.iter (fun q ->
+             Array.iter (fun p ->
+               m_pqr.{p,q,r} <- ERI.get_phys ee_ints p q r s
+             ) (Util.array_range 1 n_ao)
+           ) (Util.array_range 1 n_ao)
+         ) (Util.array_range 1 n_ao);
+      
+       (* Conversion d'un tableau a 3 indices en une matrice nao x nao^2 *)
+       let m_p_qr =
+         Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]
+         |> Bigarray.array2_of_genarray
+       in
+      
+       let m_qr_i =
+         (* (qr,i) = <i r|q s> = \sum_p <p r | q s> C_{pi} *)
+         gemm ~transa:`T ~c:m_qr_i m_p_qr m_C
+       in
+      
+       let m_q_ri =
+         (* Transformation de la matrice (qr,i) en (q,ri) *)
+         Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo)
+       in
+      
+       let m_ri_j =
+         (* (ri,j) = <i r | j s> = \sum_q <i r | q s> C_{bj} *)
+         gemm ~transa:`T ~c:m_ri_j m_q_ri m_C
+       in
+      
+       let m_r_ij =
+         (* Transformation de la matrice (ri,j) en (r,ij) *)
+         Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo)
+       in
+      
+       let m_ij_k =
+         (* (ij,k) = <i k | j s> = \sum_r <i r | j s> C_{rk} *)
+         gemm ~transa:`T ~c:m_ij_k m_r_ij m_C
+       in
+      
+       let m_ijk  =
+         (* Transformation de la matrice (ei,j) en (e,ij) *)
+         Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |]
+         |> Bigarray.array3_of_genarray
+       in
+      
+       Array.iter (fun j ->
+           Array.iter (fun i ->
+             m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});
+             m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j}  -.
+             0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_C.{s,i} +.
+                       (m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_C.{s,j})
+                 ) (Util.array_range 1 n_mo);
+                  v_d.{j} <- v_d.{j} +.  m_ijk.{j,j,j} *. m_C.{s,j}
+         ) (Util.array_range 1 n_mo)
+    ) (Util.array_range 1 n_ao);
+   
+    (Mat.init_cols n_mo n_mo ( fun i j ->
+    if i= j 
+        then 0.
+    else if  asin(  m_b12.{i,j} /. sqrt(m_b12.{i,j}**2. +. m_a12.{i,j}**2.)) > 0.
+        then 0.25 *. acos( -. m_a12.{i,j} /. sqrt(m_b12.{i,j}**2. +. m_a12.{i,j}**2.))
+        else -. 0.25 *. acos( -. m_a12.{i,j} /. sqrt(m_b12.{i,j}**2. +. m_a12.{i,j}**2.)))
+    ,Vec.sum v_d)
+  
+  in
+
+  (* 3. Function for the Boys localization like the Edminston-Ruedenberg localization *)
+  (* Coefficient matrix n_ao x n_mo -> n_mo x n_mo matrix with angle between orbitals, localization Boys like ER criterion *)
+  (* Matrix ->  matrix, float *)
+  let f_alpha_boys_er m_C = 
+      
+      let multipoles = AOBasis.multipole ao_basis in
+      let n_mo = Mat.dim2 m_C in
+      let phi_x_phi =
+          Multipole.matrix_x multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+      in    
+      let phi_y_phi =
+          Multipole.matrix_y multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+      in
+      let phi_z_phi =
+          Multipole.matrix_z multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+      in      
+      let m_b12= 
+          let b12 g =  Mat.init_cols n_mo n_mo ( fun i j ->
+                                            (g.{i,i} -. g.{j,j}) *. g.{i,j})
+          in                                
+          Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi))
+      in 
+      let m_a12 =
+          let a12 g  = Mat.init_cols n_mo n_mo ( fun i j -> 
+                        g.{i,j} *. g.{i,j} -. 0.25 *. ((g.{i,i} -. g.{j,j}) *. (g.{i,i} -. g.{j,j})))
+          in
+          Mat.add (a12 phi_x_phi) ( Mat.add (a12 phi_y_phi) (a12 phi_z_phi))
+      
+      in
+      (Mat.init_cols n_mo n_mo ( fun i j -> 
+          if i=j 
+              then 0.
+          else if  +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} >= 0.
+              then pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j}
+              else -. pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j})
+        ,0.5 *. Vec.sum(Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)))
+  in
+(* 4. Function for the original Boys localization *)
+(* Coefficient matrix n_ao x n_mo -> (n_mo x n_mo matrix with angle between orbitals, localization criterion *)
+(* Matrix -> matrix, float *)  
+  let f_theta_boys m_C = 
+      let multipoles = AOBasis.multipole ao_basis in
+      let n_mo = Mat.dim2 m_C
+      in
+      let phi_x_phi =
+          Multipole.matrix_x multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+      in    
+      let phi_y_phi =
+          Multipole.matrix_y multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+      in
+      let phi_z_phi =
+          Multipole.matrix_z multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+      in
+      let phi_x2_phi =
+          Multipole.matrix_x2 multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+      in    
+      let phi_y2_phi =
+            Multipole.matrix_y2 multipoles 
+            |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+      in
+      let phi_z2_phi =
+            Multipole.matrix_z2 multipoles 
+            |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+      in
+      let m_b_0 =
+          let b_0 g h = Mat.init_cols n_mo n_mo (fun i j ->
+              h.{i,i} +. h.{j,j} -. (g.{i,i}**2. +. g.{j,j}**2.))
+          in
+          Mat.add (b_0 phi_x_phi phi_x2_phi) (Mat.add (b_0 phi_y_phi phi_y2_phi) (b_0 phi_z_phi phi_z2_phi))  
+      in
+      let m_beta = 
+          let beta g = Mat.init_cols n_mo n_mo (fun i j ->
+              (g.{i,i} -.  g.{j,j})**2. -. 4. *. g.{i,j}**2.)
+          in
+          Mat.add (beta phi_x_phi) (Mat.add (beta phi_y_phi) (beta phi_z_phi))
+      in 
+      let m_gamma = 
+          let gamma g = Mat.init_cols n_mo n_mo (fun i j ->
+              4. *. g.{i,j} *. (g.{i,i} -. g.{j,j}) )
+          in
+           Mat.add (gamma phi_x_phi) (Mat.add (gamma phi_y_phi) (gamma phi_z_phi))
+      in 
+      let m_theta = Mat.init_cols n_mo n_mo (fun i j ->
+          if i = j 
+              then 0.
+              else +. 0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j})
+      in
+      let m_critere_B = Mat.init_cols n_mo n_mo (fun i j ->
+         0.5 *. (m_b_0.{i,j} +. 0.25 *. ((1. -. cos(4. *. m_theta.{i,j})) *. m_beta.{i,j} +. sin(4. *. m_theta.{i,j}) *. m_gamma.{i,j})))
+      
+      in
+      m_theta, Vec.sum(Mat.as_vec m_critere_B)
+  in
+  
+  (* Function to compute the new angle after rotation between the 2 orbitals j1 and j2, by the Boys like ER localization *)
+  (* n_mo x n_mo matrix with angle between orbitals, index of the orbital j1, index of the orbital j2 -> n_mo x 2 matrix with the new angles between j1 and j2, new localization criterion *)
+  (* matrix, integer, integer -> matrix, float *)
+  let f2_alpha_boys_er m_C j1 j2=
+    let multipoles = AOBasis.multipole ao_basis in
+    let n_mo = Mat.dim2 m_C
+    in
+    let phi_x_phi =
+          Multipole.matrix_x multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+    in    
+    let phi_y_phi =
+          Multipole.matrix_y multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+    in
+    let phi_z_phi =
+          Multipole.matrix_z multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+    in      
+
+    let m_b12= 
+        let b12 g =  Mat.init_cols n_mo 2 ( fun i j ->
+            if j = 1 
+                then  (g.{i,i} -. g.{j1,j1}) *. g.{i,j1}
+                else  (g.{i,i} -. g.{j2,j2}) *. g.{i,j2})
+
+        in                                
+        Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi))
+    in 
+    let m_a12 =
+        let a12 g  = Mat.init_cols n_mo 2 ( fun i j ->
+            if j = 1
+                then  g.{i,j1} *. g.{i,j1} -. 0.25 *. ((g.{i,i} -. g.{j1,j1}) *. (g.{i,i} -. g.{j1,j1}))
+                else  g.{i,j2} *. g.{i,j2} -. 0.25 *. ((g.{i,i} -. g.{j2,j2}) *. (g.{i,i} -. g.{j2,j2})))
+        in
+        Mat.add (a12 phi_x_phi) ( Mat.add (a12 phi_y_phi) (a12 phi_z_phi))
+    in
+    (Mat.init_cols n_mo 2 ( fun i j -> 
+        if i=j1 && j=1 
+            then 0.
+        else if i = j2 && j=2
+            then 0.
+        else if  +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} >= 0.
+            then pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j}
+            else -. pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j}),
+            Vec.sum(Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)))          
+  in
+  
+  (* Function to compute the new angle after rotation between the 2 orbitals j1 and j2, by the Boys localization *)
+  (* n_mo x n_mo matrix with angles between orbitals, index of the orbital j1, index of the orbital j2 -> n_mo x 2 matrix with the new angles between j1 and j2, new localization criterion *)
+  (* matrix, integer, integer -> matrix, float *)
+  let f2_theta_boys m_C j1 j2 = 
+        let multipoles = AOBasis.multipole ao_basis in
+        let n_mo = Mat.dim2 m_C
+        in
+        let phi_x_phi =
+              Multipole.matrix_x multipoles 
+              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+        in    
+        let phi_y_phi =
+              Multipole.matrix_y multipoles 
+              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+        in
+        let phi_z_phi =
+              Multipole.matrix_z multipoles 
+              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+        in
+        let phi_x2_phi =
+              Multipole.matrix_x2 multipoles 
+              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+        in    
+        let phi_y2_phi =
+              Multipole.matrix_y2 multipoles 
+              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+        in
+        let phi_z2_phi =
+              Multipole.matrix_z2 multipoles 
+              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
+        in
+        
+        let m_beta = 
+            let beta g = Mat.init_cols n_mo 2 (fun i j ->
+                if j = 1
+                    then (g.{i,i} -. g.{j1,j1})**2. -. 4. *. g.{i,j1}**2.
+                    else (g.{i,i} -. g.{j2,j2})**2. -. 4. *. g.{i,j2}**2.)
+            in Mat.add (beta phi_x_phi) (Mat.add (beta phi_y_phi) (beta phi_z_phi))
+        in 
+
+        let m_gamma = 
+            let gamma g = Mat.init_cols n_mo 2 (fun i j ->
+                if j = 1
+                   then 4. *. g.{i,j1} *. (g.{i,i} -. g.{j1,j1})
+                   else 4. *. g.{i,j2} *. (g.{i,i} -. g.{j2,j2}))
+            in Mat.add (gamma phi_x_phi) (Mat.add (gamma phi_y_phi) (gamma phi_z_phi))
+        in 
+
+        let m_theta = Mat.init_cols n_mo 2 (fun i j ->
+            if i=j1 && j=1 
+                then 0.
+            else if i = j2 && j=2
+                then 0.
+                else +. 0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j})
+        in 
+        let m_critere_B = Vec.init n_mo (fun i ->
+            phi_x2_phi.{i,i} -. phi_x_phi.{i,i}**2. +. phi_y2_phi.{i,i} -. phi_y_phi.{i,i}**2. +. phi_z2_phi.{i,i} -. phi_z_phi.{i,i}**2.) 
+        
+        in m_theta, Vec.sum( m_critere_B)
+   in
+
+  (* 5. Pattern matching : choice of the localization method *)
+  (* The method, coefficient matrix n_ao x n_mo -> n_mo x n_mo matrix with the angle between orbitals using previous function, localization criterion *) 
+  (* string, matrix -> matrix, float *) 
+
+  let m_alpha_d methode m_C = 
+    
+      let alpha methode =
+          match methode with 
+              | "Boys_er"
+              | "BOYS_ER"
+              | "boys_er" -> let alpha_boys , d_boys  = f_alpha_boys_er m_C in {m_alpha = alpha_boys; d = d_boys}
+              | "Boys"
+              | "boys"
+              | "BOYS" -> let theta_boys, b_boys = f_theta_boys m_C in {m_alpha = theta_boys; d = b_boys}
+              | "ER"
+              | "er" -> let alpha_er , d_er = f_alpha_er m_C in {m_alpha = alpha_er; d = d_er}
+              | _ -> invalid_arg "Unknown method, please enter Boys or ER"
+
+      in 
+      alpha methode
+  in
+
+  (* 5. Pattern matching : choice of the localization method to compute the n_mo x 2 new angles and not the n_mo x n_mo angles, after the rotation of i and j orbitals *)
+  (* The method, coefficient matrix n_ao x n_mo, index of the orbital i , index of the orbital j -> n_mo x 2 matrix with the new angle between orbitals i and j using previous function, new localization criterion *)  
+  (* string, matrix, integer, integer -> matrix, float *)
+  let f_speed_alpha methode m_C indice_i indice_j =
+    let alpha_speed methode =
+
+        match methode with
+            | "Boys_er"
+            | "boys_er" -> let m_alpha, critere_D = f2_alpha_boys_er m_C indice_i indice_j in {m_alpha = m_alpha; d = critere_D}
+            | "BOYS"    -> let m_alpha, critere_D = f2_theta_boys m_C indice_i indice_j in {m_alpha = m_alpha; d = critere_D}
+            | _ -> invalid_arg "Unknown method, please enter Boys or ER"
+
+      in
+      alpha_speed methode
+
+  in
+  (* 6. Norm of a matrix *)
+  (* Matrix -> float *)
+  let norm m = 
+      
+      let vec_m = Mat.as_vec m in
+      let vec2 = Vec.sqr vec_m
+      in
+      sqrt(Vec.sum vec2)
+  in 
+
+  (* 7. Calculation of the max angle[-pi/.4 , pi/.4]  and its indexes i and j in the matrix *)
+  (* n_mo x n_mo matrix with angles between orbitals, coefficient matrix n_ao x n_mo, max number of iteration -> max angle, line location, column location *)
+  (* matrix, matrix, integer -> float, integer, integer *)
+  let rec new_m_alpha m_alpha m_C n_rec_alpha=
+
+      let n_mo = Mat.dim2 m_C in
+      let alpha_m =
+          if n_rec_alpha == 0 
+              then m_alpha 
+              else 
+                  Mat.init_cols n_mo n_mo (fun i j -> 
+                  if  (m_alpha.{i,j}) >= (pi /. 2.) 
+                      then (m_alpha.{i,j} -. ( pi /. 2.))
+                  else if m_alpha.{i,j} <= -. pi /. 2.
+                      then (m_alpha.{i,j} +. ( pi /. 2.))
+                      else m_alpha.{i,j} )
+      in 
+        
+      (* Location of the max angle in the matrix *)
+      (* angle matrix -> location of the max element *)
+      (* matrix -> integer *)
+      let max_element alpha_m = 
+          Mat.as_vec alpha_m
+          |> iamax
+      in
+     
+      (* Indexes i and j of the max angle *)
+      (* -> integer, integer *)
+      let indice_ii, indice_jj =
+          let max = max_element alpha_m 
+          in
+          (max - 1) mod n_mo +1, (max - 1) / n_mo +1
+      in
+    
+      (* Value of the max angle *)
+      (* angle matrix -> value of the max angle*)
+      (* matrix -> float *)
+      let alpha alpha_m  = 
+          let i = indice_ii in
+          let j = indice_jj
+          in
+          alpha_m.{i,j}
+      in
+      let alpha_max = alpha alpha_m 
+     
+      in
+      if alpha_max < pi /. 2.
+          then {alpha_max; indice_ii; indice_jj}
+          else new_m_alpha alpha_m m_C (n_rec_alpha-1)
+  in
+
+  (* 8. Rotation matrix 2x2 *)
+  (* angle -> 2x2 rotation matrix *)
+  (* float -> matrix *)
+  let f_R alpha  =
+
+      Mat.init_cols 2 2 (fun i j -> 
+          if i=j 
+              then cos alpha
+          else if i>j 
+              then  sin alpha 
+              else -. sin alpha )
+  in
+
+  (* 9. Rotation matrix nxn *)
+  (* Coefficient matrix n_ao x n_mo, angle -> new coefficient matrix n_ao x n_mo after a rotation with a small angle *)
+  (* matrix, float -> matrix *)
+  let rotate m_C angle=
+ 
+      let n_mo =  Mat.dim2 m_C in
+
+      (* Antisymmetrization of the matrix and diagonal terms = 0, X *)
+      (* angle -> Antisymmetric matrix with this angle and diagonal terms = 0 *)
+      (* float -> matrix *)
+      let m_angle angle=
+
+          Mat.init_cols n_mo n_mo (fun i j ->
+              if i = j 
+                  then 0.
+              else if i>j 
+                  then -. angle
+                  else angle)
+      in
+  
+      let m = m_angle angle in
+  
+      (* Square of the matrix X² *)
+      (* -> X . X *)
+      (* -> matrix*)
+      let mm = gemm m m in
+
+      (* Diagonalization of X² *)
+      (* diag = vector that contains the eigenvalues (-tau²), mm contains the eigenvectors (W) *)
+      (* -> vector  *)
+      let diag = syev ~vectors:true mm in
+  
+      (* -tau² to tau² *)
+      (* -> vector that contains the eigenvalues (tau²) *)
+      (* -> vector *)
+      let square_tau= Vec.abs diag in
+  
+      (* tau² to tau *)
+      (* -> vector that contains the eigenvalues tau *)
+      (* -> vector *)
+      let tau = Vec.sqrt square_tau in 
+  
+      (* Calculation of cos tau from the vector tau *)
+      (* -> cos(tau) *)
+      (* -> integer *)
+      let cos_tau = Vec.cos tau in
+  
+      (* Matrix cos tau *)
+      (* -> matrix with diagonal terms cos(tau) and 0 also *)
+      (* -> matrix *)
+      let m_cos_tau =
+ 
+          Mat.init_cols n_mo n_mo (fun i j ->
+              if i=j
+                  then cos_tau.{i}
+                  else 0.)
+      in
+
+      (* Calcul of sin(tau) from the vector tau *)
+      let sin_tau = Vec.sin tau in 
+      
+      (* Matrix sin tau *)
+      (* -> matrix with diagonal terms sin(tau) and 0 also *)
+      (* -> matrix *)
+      let m_sin_tau =
+
+          Mat.init_cols n_mo n_mo (fun i j ->
+              if i=j 
+                  then sin_tau.{i} 
+                  else 0.) 
+      in
+
+      (* Calculation of the transposed matrix of W (X²) *)
+      (* -> transposed matrix of W *)
+      (* -> matrix *)
+      let transpose_mm = Mat.transpose_copy mm in
+
+      (* Calculation of the vector tau^-1 *)
+      (* -> vector tau^-1 *)
+      (* -> vector *)
+      let tau_1 =
+
+          Vec.init n_mo (fun i ->
+              if tau.{i}<= 0.001
+                  then 1.
+                  else 1. /. tau.{i})
+      in
+      (* Calculation of the matrix tau^⁻1 *)
+      (* -> matrix tau^-1 *)
+      (* -> matrix *)
+      let m_tau_1 = 
+
+          Mat.init_cols n_mo n_mo (fun i j -> 
+              if i=j 
+                  then tau_1.{i} 
+                  else 0.) 
+      in
+
+      (* gemm mm (gemm m_cos_tau transpose_mm) -> W cos(tau) Wt *)
+      (* -> matrix *)
+      let a = gemm mm (gemm m_cos_tau transpose_mm) in
+
+      (* gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) -> W tau^-1 sin(tau) Wt X *)
+      (* -> matrix *)
+      let b = gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) in
+
+      (* Sum of a + b -> Rotation matrix *)
+      (* -> Rotation matrix *)
+      (* -> matrix *)
+      let m_r = Mat.add a b 
+      in
+      gemm m_C m_r
+  in
+
+  (* Function to extract 2 vectors i and j in a matrix *)
+  (* Coefficient matrix n_ao x n_mo, index of the orbital i,  index of the orbital j -> n_mo x 2 matrix with the eigenvetors i and j*)
+  (* matrix, integer, integer -> matrix *)
+  let f_Ksi indice_i indice_j m_C  =
+      let n_ao = Mat.dim1 m_C 
+
+      in
+      Mat.init_cols n_ao 2 (fun i j -> 
+          if j=1 
+              then m_C.{i,indice_i} 
+              else m_C.{i,indice_j} )
+  in
+
+  (* Function to apply a rotation with a 2x2 rotation matrix on a n_mo x 2 matrix that contains i and j *)
+  (* 2x2 rotation matrix, n_mo x 2 eigenvetors(i and j) matrix -> n_mo x 2 new eigenvetors(i~ and j~) after rotation *)
+  (* matrix, matrix -> matrix *)
+  let f_Ksi_tilde m_R m_Ksi = gemm m_Ksi m_R 
+  in
+
+  (* Function to create intermediate coefficient matrix in order to delete the i j orbitals and add i~ j~ in the coefficient matrix *)
+  (* n_mo x 2 matrix with k and l orbitals, index of k, index of l, n_ao x n_mo matrix coefficient -> matrix where all terms are equal to 0 except the k and l columns that contain the eigenvetors k and l *)
+  (* matrix, integer, integer, matrix -> matrix *)
+  let f_k mat indice_i indice_j m_C =
+      let n_mo = Mat.dim2 m_C in
+      let n_ao = Mat.dim1 m_C
+      
+      in
+      Mat.init_cols n_ao n_mo (fun i j -> 
+          if j=indice_i 
+              then mat.{i,1}
+          else if j=indice_j 
+              then mat.{i,2}
+          else 0.)
+  in
+
+  (* Function to create a intermediate angle matrix in order to delete all the previous angle between i and all the other orbitals, j and all the other orbitals.
+  And add the new angles between these orbitals *)
+  (* n_mo x 2 matrix with angles between k and all the other orbitals l and all the other orbitals, index of k, index of l, coefficent matrix n_ao x n_mo -> matrix where all terms are equal to 0 except the k and l lines and columns that contains the angles between the orbitals *)
+  (* matrix, integer, integer, matrix -> matrix  *)
+  let f_k_angle mat indice_i indice_j m_C = 
+      let n_mo = Mat.dim2 m_C
+      
+      in
+      Mat.init_cols n_mo n_mo (fun i j -> 
+          if j=indice_i 
+              then mat.{i,1}
+          else if j=indice_j 
+              then mat.{i,2}
+          else if i = indice_i
+              then -. mat.{j,1}
+          else if i = indice_j
+              then -. mat.{j,2}
+          else 0.)
+  in
+
+  (* Intermediate matrix where the i and j orbitals are equal to 0 *)
+  (* Coefficient matrix n_ao x n_mo, matrix where terms are equal to 0 except the k and l columns that contain the eigenvetors k and l *) 
+  (* matrix, matrix -> matrix *)
+  let f_interm m_C mat = 
+
+      Mat.sub m_C mat
+  in
+
+  (* Function to compute the new coefficient matrix after a rotation between 2 orbitals *)
+  (* coefficient matrix n_ao x n_mo, angle matrix, parameter -> new coefficient matrix n_ao x n_mo *)
+  let new_m_C m_C m_alpha epsilon =
+
+      (* alphaij contains the max angle of the angle matrix and the indexes *)
+      let n_rec_alpha = 10 in 
+      let alphaij = new_m_alpha m_alpha m_C n_rec_alpha in
+
+      (* Value of the angle *)
+      let alpha = (alphaij.alpha_max) *. epsilon in
+
+      (* Localtion of the max angle *)
+      let indice_i = alphaij.indice_ii in
+      let indice_j = alphaij.indice_jj in
+
+      (* Rotation matrix *)
+      let m_R = f_R alpha in
+
+      (* Matrix with i and j orbitals *)
+      let m_Ksi = f_Ksi indice_i indice_j m_C in
+
+      (* Matrix with the new i~ and j~ orbitals *)
+      let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi in
+
+      (* Matrix to delete the i and j orbitals in the coefficient matrix *)
+      let m_Psi = f_k m_Ksi indice_i indice_j m_C in
+
+      (* Matrix to add the i~ and j~ orbitals in the coefficient matrix *)
+      let m_Psi_tilde = f_k m_Ksi_tilde indice_i indice_j m_C in
+
+      (* Coefficient matrix without the i and j orbitals *)
+      let m_interm = f_interm m_C m_Psi
+
+      in
+      (* Matrice after rotation, max angle, index i, index j *)
+      (Mat.add m_Psi_tilde m_interm, alpha, indice_i, indice_j)
+  in
+
+  (* Function to compute the new angle matrix after one rotation between 2 orbitals *)
+  (* Previous angle matrix n_mo x n_mo, new coefficient matrix n_ao x n_mo, previous angle matrix, index of the orbital i, index of the orbital j -> new angle matrix, new localization criterion *)
+  (* matrix, matrix, matrix, integer, integer -> matrix, float *)
+  let m_alpha m_new_m_C prev_m_alpha indice_i indice_j =
+      
+      (* n_mo x 2 matrix that contains the new angles between i,j with the other orbitals *)    
+      let speed_alphad = f_speed_alpha methode m_new_m_C indice_i indice_j in
+
+      (* New localization criterion *)
+      let critere_D = speed_alphad.d in
+
+      (* Previous angles *)
+      let ksi_angle = f_Ksi indice_i indice_j prev_m_alpha in
+      
+      (* New angles *)
+      let ksi_angle_tilde = speed_alphad.m_alpha in
+
+      (* Matrix to delete the previous angles *)
+      let psi_angle = f_k_angle ksi_angle indice_i indice_j m_new_m_C in
+
+      (* Matrix to add the new angles *)
+      let psi_tilde_angle = f_k_angle ksi_angle_tilde indice_i indice_j m_new_m_C in
+
+      (* Matrix without the angles between i, j and the orbitals *)
+      let m_interm = f_interm prev_m_alpha psi_angle
+
+      in
+      (* New angle matrix, new localization criterion *)
+      (Mat.add m_interm psi_tilde_angle, critere_D)
+  in
+
+
+
+(* III. Localization *)
+
+  (* Loop to compute the new coefficient matrix avec n orbital rotations *)
+  (* Coefficeient matrix n_ao x n_mo, localization method, parameter, max number of iteration, previous localization criterion, previous angle matrix, convergence criterion -> new coefficient matrix after n orbital rotations *)
+  (* matrix, string, float, integer, float, matrix, float -> matrix *) 
+  let rec localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc =
+
+      if n == 0 
+        then m_C
+        else
+
+           (* New coefficient matrix after one rotation *)
+           let m_new_m_C, alpha_max, indice_i, indice_j = new_m_C m_C prev_m_alpha epsilon in
+          
+           (* New angle matrix after one rotation *)
+           let m_alpha, critere_D = m_alpha m_new_m_C prev_m_alpha indice_i indice_j  in
+
+           (* Matrix norm, angle average  *)
+    
+           let norm_alpha = norm m_alpha in
+           let moyenne_alpha = norm_alpha /. (float_of_int((Mat.dim1 m_C)* (Mat.dim2 m_C))) in
+           let alpha = alpha_max /. epsilon in
+    
+           Printf.printf "%i %f %f %f %f\n%!" (iteration-n) critere_D alpha norm_alpha moyenne_alpha;
+   
+           (* Convergence criterion *)
+           if alpha_max**2.  < cc**2.
+               then m_new_m_C
+               else localisation m_new_m_C methode epsilon (n-1) critere_D m_alpha cc
+  in
+
+  (* 2. Localization function *)
+  let localise localisation m_C methode epsilon n cc =
+
+      let alphad = m_alpha_d methode m_C in
+      let prev_m_alpha = alphad.m_alpha in 
+      let prev_critere_D = 0. 
+      in
+      localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc
+  in
+
+(* IV : Functions to split and assemble matrix *)
+
+  (* Function to create a integer list from a vector of float *)
+  (* float vector -> integer list *)
+  let int_list vec = 
+
+      let float_list = Vec.to_list vec in
+      let g a = int_of_float(a)
+      in
+      List.map g float_list
+  in
+
+  (* Function to create a list from the missing elements of an other list*)
+  (* Coefficient matrix n_ao x n_mo, list -> list of missing element of the previous list *)
+  (* matrix, integer list -> integer list *)
+  let miss_elem mat list = 
+
+      let n_mo = Mat.dim2 mat in
+      let vec = Vec.init (n_mo) (fun i ->
+          if List.mem i list
+              then 0.
+              else float_of_int(i))
+      in
+      let list_int = int_list vec    
+      in
+      List.filter (fun x -> x > 0) list_int
+  in
+
+  (* Function to split a matrix in 2 matrix, the first matrix corresponds to the column number whose number is in the list, the second matrix corresponds to the column which are not in the list *)
+  (* Coefficient matrix n_ao x n_mo, list -> matrix, matrix *)
+  (* matrix, integer list -> matrix, matrix*)
+  let split_mat mat list =
+
+      let vec_of_mat = Mat.to_col_vecs mat in
+      let f a = vec_of_mat.(a-1) in
+      let vec_list_1 = List.map f list in
+      let list_2 = miss_elem mat list in
+      let vec_list_2 = List.map f list_2
+      in
+      (Mat.of_col_vecs_list vec_list_1,Mat.of_col_vecs_list vec_list_2)
+  in
+
+
+  (* Function to create a matrix from a list*)
+  (* Coefficient matrix n_ao x n_mo, list of orbitals -> matrix with these orbitals *)
+  (* matrix, integer list -> matrix*) 
+  let create_mat mat list = 
+
+      let vec_of_mat = Mat.to_col_vecs mat in
+      let f a = vec_of_mat.(a-1) in
+      let vec_list = List.map f list 
+      in
+      Mat.of_col_vecs_list vec_list
+  in
+
+  (* List of the occupied MOs *)
+  let list_occ = 
+
+      let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))
+      in
+      int_list vec_occ
+  in
+
+  (* List of the core MOs *)
+  let list_core = 
+      let vec_core = Vec.init (n_core) (fun i -> float_of_int(i))
+      in int_list vec_core
+  in
+  
+  (* Function to create a list of pi MOs *)
+  (* Coefficient matrix n_ao x n_mo -> list of pi orbitals *)
+  (* matrix -> integer list *)
+  let list_pi mat =
+      let n_mo = Mat.dim2 mat in
+      let m_pi = Mat.init_cols n_ao (n_mo+1) ( fun i j ->
+          if j=1 && m_C.{i,j}**2. < 10e-8**2.
+              then 0.
+          else if j=1 
+              then m_C.{i,j}
+          else if mat.{i,j-1}**2. < 10e-8**2.
+              then 0.
+              else mat.{i,j-1})
+      in
+      let vec = Mat.to_col_vecs m_pi in
+      let vec_dot = Vec.init ((Mat.dim2 mat)) (fun i -> dot vec.(0) vec.(i)) in
+      let vec_pi = Vec.init ((Mat.dim2 mat)) (fun i ->
+          if vec_dot.{i} = 0.
+              then float_of_int(i)
+              else 0.)
+      in
+      let list_pi = int_list vec_pi in
+      let rec remove x list =
+          match list with
+            | [] -> []
+            | var :: tail -> if var = x
+                       then remove x tail
+                       else var :: (remove x tail)
+      in
+      remove 0 list_pi
+  in 
+
+  (* Function to create a list of the (n-x_end) first MOs of a matrix *)
+  (* Coefficient matrix n_ao x n_mo -> list of the (n-x_end) first MOs of the matrix *)
+  (* matrix -> integer list *)
+  let list_x_end mat =
+      let n = Mat.dim2 mat in
+      let vec_x_end = Vec.init (n-x_end)  (fun i -> float_of_int(i))
+      in
+      int_list vec_x_end
+  in
+
+  (* Function to create a list of the (n/2) first MOs of a matrix *)
+  (* Coefficient matrix n_ao x n_mo -> list of the (n/2) first MOs of the matrix *)
+  (* matrix -> integer list *)
+  let list_s_12 mat =
+      let n = Mat.dim2 mat in
+      let vec_12 = Vec.init (n/2)  (fun i -> float_of_int(i))
+      in
+      int_list vec_12
+  in
+   
+  (* Function to assemble to matrix with the same number of line *)
+  (* Coefficient matrix n_ao x n, coefficient matrix n_ao x m -> Coefficient matrix n_ao x (n+m) *)
+  (* matrix, matrix -> matrix*) 
+  let assemble mat_1 mat_2 =
+
+      let v_mat_1 = Mat.to_col_vecs mat_1 in
+      let v_mat_2 = Mat.to_col_vecs mat_2 in
+      let m = Array.append v_mat_1 v_mat_2 
+      in
+      Mat.of_col_vecs m
+  in
+
+(* V. Calculation  *)
+
+   (* 1. Caption *)
+  let text = "Method : " ^ methode ^ "; Separation type of the valence occupied MOs : " ^ type_separation_occu  ^ "; Separation type of the virtual MOs : " ^ type_separation_vir ^ "; Max number of iteration : " ^ string_of_int(iteration) ^ "; cc : " ^ string_of_float(cc)^"; epsilon : "^string_of_float(epsilon)^"; charge : "^string_of_int(charge)^"; Rotation pre angle : "^string_of_float(pre_angle)
+  in
+  let caption = "n Criterion max_angle norm_angle average_angle"
+  in
+  Printf.printf "%s\n" "";
+  Printf.printf "%s\n" "Orbitals localization ";
+  Printf.printf "%s\n" "";
+  Printf.printf "%s\n" text;
+  Printf.printf "%s\n" "";
+  Printf.printf "%s\n" caption;
+  Printf.printf "%s\n" "";
+
+  (* Separation of the occupied and virtual MOs, valence and core MOs *)
+  let m_occ , m_vir = split_mat m_C list_occ in
+   
+  let m_core, m_occu = split_mat m_occ list_core in
+
+  (* Localization function *) 
+  (* Coefficient matrix n_ao x n_mo -> localized coefficient matrix n_ao x n_mo *)
+  (* matrix -> matrix *)
+  let f_localise mat =
+
+      if Mat.dim2 mat = 0
+          then Printf.printf "%s\n" "0 orbital, no localization";
+
+      let new_mat =
+          if Mat.dim2 mat = 0
+              then mat
+              else rotate mat pre_angle
+      in
+      if Mat.dim2 mat = 0
+          then new_mat 
+          else localise localisation new_mat methode epsilon iteration cc
+  in
+  
+  (* Core MOs localization *)
+  let loc_m_core = 
+          Printf.printf "%s\n" "I/III"; 
+          Printf.printf "%s\n" "Localization of core orbitals";
+          Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 m_core);
+          Printf.printf "%s\n" "";
+          
+          f_localise m_core in
+          
+          Printf.printf "%s\n" "";
+          Printf.printf "%s\n" "End I/III";
+          Printf.printf "%s\n" "";
+
+  (* Localization function without separation *)
+  let f_loc_assemble_1 mat =
+      
+      Printf.printf "%s\n" "";
+      Printf.printf "%s\n" "Localization 1/1";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat);
+      
+      f_localise mat in
+      
+      Printf.printf "%s\n" "";
+
+  (* Localization function with 1 MOs separation *)
+  let f_loc_assemble_12 mat list_12  =
+      
+      if List.length(list_12 mat) + List.length(miss_elem mat (list_12 mat)) <> (Mat.dim2 mat) 
+          then Printf.printf "%s\n" "Bug : List length problem (list_12 or miss_elem list_12)";
+  
+      let mat_1, mat_2 = split_mat mat (list_12 mat) in
+
+      Printf.printf "%s\n" "";
+      Printf.printf "%s\n" "Localization 1/2";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_1);
+      let m_loc_1 = f_localise mat_1 in
+      Printf.printf "%s\n" "";
+
+      Printf.printf "%s\n" "Localization 2/2";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_2);
+      let m_loc_2 = f_localise mat_2 in
+      Printf.printf "%s\n" "";
+
+      assemble m_loc_1 m_loc_2
+  in
+
+  (* Localization function with 2 MOs separations *)
+  let f_loc_assemble_123 mat list_12 list_2ab =
+
+      if List.length(list_12 mat) + List.length(miss_elem mat (list_12 mat)) <> (Mat.dim2 mat)
+          then Printf.printf "%s\n" "Bug : List length problem (list_12 or miss_elem list_12)";
+
+      let mat_1, mat_2 = split_mat mat (list_12 mat) in
+
+      if List.length(list_2ab mat_2) + List.length(miss_elem mat_2 (list_2ab mat_2)) <> (Mat.dim2 mat_2)
+          then Printf.printf "%s\n" "Bug : List length problem (list_2ab or miss_elem list_2ab)";
+      
+      let mat_2a, mat_2b = split_mat mat_2 (list_2ab mat_2) in
+ 
+      Printf.printf "%s\n" "";
+      Printf.printf "%s\n" "Localization 1/3";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_1);
+      let m_loc_1 = f_localise mat_1 in
+      Printf.printf "%s\n" "";
+
+      Printf.printf "%s\n" "Localization 2/3";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_2a);
+      let m_loc_2a = f_localise mat_2a in
+      Printf.printf "%s\n" "";
+
+      Printf.printf "%s\n" "Localization 3/3";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_2b);
+      let m_loc_2b = f_localise mat_2b in
+      Printf.printf "%s\n" "";
+
+      let m_loc_2ab = assemble m_loc_2a m_loc_2b
+      in
+      
+      assemble m_loc_1 m_loc_2ab
+  in
+
+  (* Localization function with 3 MOs separations *)
+  let f_loc_assemble_1234 mat list_12 list_1ab list_2ab  =
+
+      if List.length(list_12 mat) + List.length(miss_elem mat (list_12 mat)) <> (Mat.dim2 mat)
+          then Printf.printf "%s\n" "Bug : List length problem (list_12 or miss_elem list_12)";
+
+      let mat_1, mat_2 = split_mat mat (list_12 mat) in
+
+      if List.length(list_1ab mat_1) + List.length(miss_elem mat (list_1ab mat_1)) <> (Mat.dim2 mat_1)
+          then Printf.printf "%s\n" "Bug : List length problem (list_1ab or miss_elem list_1ab)";
+
+      let mat_1a, mat_1b = split_mat mat_1 (list_1ab mat_1) in
+
+      if List.length(list_2ab mat_2) + List.length(miss_elem mat (list_2ab mat_1)) <> (Mat.dim2 mat_2)
+          then Printf.printf "%s\n" "Bug : List length problem (list_2ab or miss_elem list_2ab)";
+
+      let mat_2a, mat_2b = split_mat mat_2 (list_2ab mat_2) in
+
+      Printf.printf "%s\n" "";
+      Printf.printf "%s\n" "Localization 1/4";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_1a);
+      let m_loc_1a = f_localise mat_1a in
+      Printf.printf "%s\n" "";
+
+      Printf.printf "%s\n" "Localization 2/4";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_1b);
+      let m_loc_1b = f_localise mat_1b in
+      Printf.printf "%s\n" "";
+
+      Printf.printf "%s\n" "Localization 3/4";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_2a);
+      let m_loc_2a = f_localise mat_2a in
+      Printf.printf "%s\n" "";
+
+      Printf.printf "%s\n" "Localization 4/4";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_2b);
+      let m_loc_2b = f_localise mat_2b in
+      Printf.printf "%s\n" "";
+
+      let m_loc_1ab = assemble m_loc_1a m_loc_1b in
+      let m_loc_2ab = assemble m_loc_2a m_loc_2b
+      in
+      assemble m_loc_1ab m_loc_2ab
+  in
+
+  (* Localization function with 4 separations "by hands" *)
+  let f_loc_assemble_main mat list_1 list_2 list_3 list_4 =
+      
+      if (List.length(list_1) + List.length(list_2) + List.length(list_3) + List.length(list_4)) <> (Mat.dim2 mat)
+          then Printf.printf "%s\n" "Bug : List length problem";
+
+      let mat_1 = create_mat mat list_1 in
+      let mat_2 = create_mat mat list_2 in
+      let mat_3 = create_mat mat list_3 in
+      let mat_4 = create_mat mat list_4 in
+      
+      Printf.printf "%s\n" "";
+      Printf.printf "%s\n" "Localization 1/4";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_1);
+      let loc_mat_1 = f_localise mat_1 in
+      Printf.printf "%s\n" "";
+
+      Printf.printf "%s\n" "Localization 2/4";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_2);
+      let loc_mat_2 = f_localise mat_2 in
+      Printf.printf "%s\n" "";
+
+      Printf.printf "%s\n" "Localization 3/4";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_3);
+      let loc_mat_3 = f_localise mat_3 in
+      Printf.printf "%s\n" "";
+
+      Printf.printf "%s\n" "Localization 4/4";
+      Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_4);
+      let loc_mat_4 = f_localise mat_4 
+      
+      in
+      Printf.printf "%s\n" "";
+      assemble (assemble loc_mat_1 loc_mat_2) (assemble loc_mat_3 loc_mat_4)
+  in
+  
+  (* Pattern matching for the separation type of the occupied MOs *)
+  let m_assemble_occ, m_assemble_occu =
+
+      let m_assemble_occu =
+          Printf.printf "%s\n" "II/III";
+          Printf.printf "%s\n" "Localization of valence orbitals";
+          Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 m_occu);
+          match type_separation_occu with
+              | "full"         -> let m_assemble_occu = f_loc_assemble_1 m_occu in m_assemble_occu
+              | "pi/sigma"     -> let m_assemble_occu = f_loc_assemble_12 m_occu list_pi in m_assemble_occu
+              | "pi/sigma/end" -> let m_assemble_occu = f_loc_assemble_123 m_occu list_pi list_x_end in m_assemble_occu
+              | "pi/1/2"       -> let m_assemble_occu = f_loc_assemble_123 m_occu list_pi list_s_12 in m_assemble_occu
+              | "1/2"          -> let m_assemble_occu = f_loc_assemble_12 m_occu list_s_12 in m_assemble_occu 
+              | "pi1/pi2/1/2"  -> let m_assemble_occu = f_loc_assemble_1234 m_occu list_pi list_s_12 list_s_12 in m_assemble_occu
+              | "By_hands"    -> let m_assemble_occu = f_loc_assemble_main m_occu (miss_elem m_occ []) [] [] [] in m_assemble_occu
+              | _              -> invalid_arg "Unknown separation type of valence orbitals"
+      in
+      Printf.printf "%s\n" "End II/III";
+      Printf.printf "%s\n" "";
+      (assemble loc_m_core m_assemble_occu, m_assemble_occu)
+  in
+
+  (* Pattern matching for the separation type of the virtual MOs *)
+  let m_assemble_vir =
+       
+      let m_assemble_vir =
+          Printf.printf "%s\n" "III/III";
+          Printf.printf "%s\n" "Localization of virtual orbitals";
+          Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 m_vir);
+          match type_separation_vir with
+              | "full"         -> let m_assemble_vir = f_loc_assemble_1 m_vir in m_assemble_vir
+              | "pi/sigma"     -> let m_assemble_vir = f_loc_assemble_12 m_vir list_pi in m_assemble_vir
+              | "pi/sigma/end" -> let m_assemble_vir = f_loc_assemble_123 m_vir list_pi list_x_end in m_assemble_vir
+              | "pi/1/2"       -> let m_assemble_vir = f_loc_assemble_123 m_vir list_pi list_s_12 in m_assemble_vir
+              | "1/2"          -> let m_assemble_vir = f_loc_assemble_12 m_vir list_s_12 in m_assemble_vir
+              | "pi1/pi2/1/2"  -> let m_assemble_vir = f_loc_assemble_1234 m_vir list_pi list_s_12 list_s_12 in m_assemble_vir
+              | "By_hands"    -> let m_assemble_vir = f_loc_assemble_main m_vir (miss_elem m_vir []) [] [] [] in m_assemble_vir
+              | _              -> invalid_arg "Unknown separation type of virtual orbitals"
+      in
+      Printf.printf "%s\n" "End III/III";
+      Printf.printf "%s\n" "";
+      m_assemble_vir
+  in
+
+  (* Tack occupied and virtual MOs together*)
+  let m_assemble_loc = assemble m_assemble_occ m_assemble_vir
+  in
+
+  (* VI. Analysis *)
+
+  (* Function to compute the spatial extent of MOs *)
+  (* Coefficient matrix n_ao x n_mo -> vector where the i th element is equal to the spatial extent of the i th orbital *)
+  (* Matrix -> vector *)
+  let distrib mat= 
+      if Mat.dim1 mat = 0
+      then Mat.as_vec mat 
+      else
+      (let multipoles = AOBasis.multipole ao_basis in
+      let n_mo = Mat.dim2 mat in
+      let phi_x_phi =
+          Multipole.matrix_x multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
+      in    
+      let phi_y_phi =
+          Multipole.matrix_y multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
+      in
+      let phi_z_phi =
+          Multipole.matrix_z multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
+      in
+      let phi_x2_phi =
+          Multipole.matrix_x2 multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
+      in    
+      let phi_y2_phi =
+          Multipole.matrix_y2 multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
+      in
+      let phi_z2_phi =
+          Multipole.matrix_z2 multipoles 
+          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
+      in 
+      Vec.init n_mo (fun i ->
+      phi_x2_phi.{i,i} -. phi_x_phi.{i,i}**2. +. phi_y2_phi.{i,i} -. phi_y_phi.{i,i}**2. +. phi_z2_phi.{i,i} -. phi_z_phi.{i,i}**2.))
+  in
+ 
+  (* Sorting function *) 
+  (* tri insertion (<) list -> list in ascending order *)
+  (* parameter, list -> list *)
+  let rec insere comparaison elem liste = match liste with
+          |  [] -> elem::[]
+          |  tete::queue -> if comparaison elem tete
+                                then elem :: liste
+                                else tete :: insere comparaison elem queue
+  in
+  let rec tri_insertion comp = function
+          |  [] -> []
+          |  tete::queue -> insere comp tete (tri_insertion comp queue)
+  in
+
+  (* Sorting function to sort orbital spatial extent by ascending order *)
+  (* coefficient matrix n_ao x n_mo -> matrix n_ao x 2, first column corresponds to the spatial extent in function of the MOs number, the second column corresponds to the spatial extent in ascending order *)
+  (* matrix -> matrix *)
+  let m_distribution mat = 
+      let vec_distrib= distrib mat in
+      let list_tri_distrib = tri_insertion (<) (Vec.to_list vec_distrib) in
+      let vec_tri_distrib = Vec.of_list list_tri_distrib
+      in
+      Mat.init_cols (Vec.dim vec_distrib) 2 (fun i j ->
+              if j=1 
+                  then vec_distrib.{i}
+                  else vec_tri_distrib.{i})
+
+  in
+
+  (* Average/variance/standart deviation/median of spatial extent *)
+  (* Function to compute the average *)
+  let f_average mat =
+          if Mat.dim1 mat = 0
+          then 0.
+          else  Vec.sum (distrib mat) /. float_of_int(Vec.dim (distrib mat))
+  in
+
+  (* Function to compute the variance *)
+  let f_variance mat =
+          if (Vec.dim (distrib mat)) = 0
+          then 0.
+          else  Vec.sum(Vec.init (Vec.dim (distrib mat)) (fun i ->
+     ((distrib mat).{i}-. (f_average mat))**2.)) /. float_of_int(Vec.dim (distrib mat))
+  in
+
+  (* Function to compute the standard deviation *)
+  let f_stand_dev mat = 
+          if (Vec.dim (distrib mat)) = 0
+          then 0.
+          else  sqrt(abs_float(f_variance mat)) 
+  in
+  
+  (* Fonction de calcul de la mediane *)
+  let f_median mat = 
+          if (Vec.dim (distrib mat)) = 0
+          then 0.
+          else
+      let vec_distrib= distrib mat in
+      let list_tri_distrib = tri_insertion (<) (Vec.to_list vec_distrib) in
+      let vec_tri_distrib = Vec.of_list list_tri_distrib
+      in
+      if (Vec.dim vec_distrib) mod 2 = 0
+          then (vec_tri_distrib.{(Vec.dim vec_distrib)/2} +. vec_tri_distrib.{((Vec.dim vec_distrib)/2)+1}) /. 2. 
+          else vec_tri_distrib.{int_of_float(float_of_int(Vec.dim vec_distrib)/. 2.) + 1}
+  in
+  (* Display function *)
+  (* Coefficient matrix n_ao x n_mo -> string *)
+  (* matrix -> string *)
+  let analyse mat = 
+      let average = f_average mat in
+      let variance = f_variance mat in
+      let stand_dev = f_stand_dev mat in
+      let median = f_median mat
+      in
+      "Average : "^string_of_float(average)^"; Median : "^string_of_float(median)^"; Variance : "^string_of_float(variance)^"; Standard deviation : "^string_of_float(stand_dev)
+  in
+  (* Display the numer of MOs *)
+  (* Coefficient matrix n_ao x n_mo -> string *)
+  (* matrix -> string *)
+  let n_orb mat = "Number of molecular orbitals : "^string_of_int(Mat.dim2 mat) in
+
+  (* Display *)
+  Printf.printf "%s\n" "";
+  Printf.printf "%s\n" (n_orb m_C);
+  Printf.printf "%s %i\n" "Number of pi orbitals : " (List.length (list_pi m_C));
+  Printf.printf "%s %i\n" "Number of occupied pi orbitals : " (List.length (list_pi m_occ));
+  Printf.printf "%s %i\n" "Number of virtual pi orbitals : " (List.length (list_pi m_vir));
+  Printf.printf "%s %i\n" "Number of sigma orbitals : " (List.length (miss_elem m_C (list_pi m_C)));
+  Printf.printf "%s %i\n" "Number of occupied sigma orbitals : " (List.length (miss_elem m_occ (list_pi m_occ)));
+  Printf.printf "%s %i\n" "Number of virtual sigma orbitals : " (List.length (miss_elem m_vir (list_pi m_vir)));
+
+  Printf.printf "%s\n" "";
+
+  Printf.printf "%s %s\n" "All MOs before the localization" (analyse m_C);
+  Printf.printf "%s %s\n" "All MOs after the localization" (analyse m_assemble_loc);
+  Printf.printf "%s\n" "";
+  Util.debug_matrix "Distribution of the spatial extent, before localization / sorting before localization / after localization / sorting after localization" (assemble(m_distribution m_C) (m_distribution m_assemble_loc));
+  
+  Printf.printf "%s\n" "";
+  Printf.printf "%s %s\n" "Occupied orbitals : " (n_orb m_occ);
+  Printf.printf "%s %s\n" "Occupied orbitals before localization"  (analyse m_occ);
+  Printf.printf "%s %s\n" "Occupied orbitals after localization" (analyse m_assemble_occ);
+  Printf.printf "%s\n" "";
+
+  Printf.printf "%s\n" "";
+  Printf.printf "%s %s\n" "Core orbitals : " (n_orb m_core);
+  Printf.printf "%s %s\n" "Core orbitalses before localization : " (analyse m_core);
+  Printf.printf "%s %s\n" "Core orbitals after localization : " (analyse loc_m_core);
+  Printf.printf "%s\n" "";
+
+  Printf.printf "%s\n" "";
+  Printf.printf "%s %s\n" "Valence orbitals : " (n_orb m_occu);
+  Printf.printf "%s %s\n" "Valence orbitals before localization : " (analyse m_occu);
+  Printf.printf "%s %s\n" "Valence rbitals after localization : " (analyse m_assemble_occu);
+  Printf.printf "%s\n" "";
+
+  Printf.printf "%s\n" "";
+  Printf.printf "%s %s\n" "Virtual orbitals : " (n_orb m_vir);
+  Printf.printf "%s %s\n" "Virtual orbitals before localization : " (analyse m_vir);
+  Printf.printf "%s %s\n" "Virtual orbitals after localization : " (analyse m_assemble_vir);
+  Printf.printf "%s\n" "";
+
+m_assemble_loc 
+
+
diff --git a/run_Localisation.ml b/run_Localisation.ml
new file mode 100644
index 0000000..32eb5ae
--- /dev/null
+++ b/run_Localisation.ml
@@ -0,0 +1,124 @@
+(* ./run_Localisation.native -b basis -x xyz -c charge -o name.mos -i EZFIOdirectory > name.out *)
+(* 
+#.mos contains the localised orbitales 
+#.out contains the localization convergence and the analysis of the spatial extent of the orbitales
+*)
+
+open Lacaml.D
+
+let read_qp_mo dirname =
+  let ic = Unix.open_process_in ("zcat "^dirname^"/mo_basis/mo_coef.gz") in
+  let check = String.trim (input_line ic) in
+  assert (check = "2");
+  let int_list = 
+    input_line ic
+    |> String.split_on_char ' ' 
+    |> List.filter (fun x -> x <> "")
+    |> List.map int_of_string
+  in
+  let n_ao, n_mo = 
+    match int_list with
+    | [ x ; y ] -> x, y
+    | _ -> assert false
+  in
+  let result = 
+          Mat.init_cols n_ao n_mo (fun i j ->
+            let s = input_line ic in
+            Scanf.sscanf s " %f " (fun x -> x)
+          )
+  in
+  let exit_code = Unix.close_process_in ic in
+  assert (exit_code = Unix.WEXITED 0);
+  result
+
+
+
+let () =
+  let open Command_line in
+  begin
+    set_header_doc (Sys.argv.(0) ^ " - QCaml command");
+    set_description_doc "Localizes MOs";
+    set_specs
+      [ { short='b' ; long="basis" ; opt=Mandatory;
+          arg=With_arg "<string>";
+          doc="Name of the file containing the basis set"; } ;
+
+        { short='x' ; long="xyz" ; opt=Mandatory;
+          arg=With_arg "<string>";
+          doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
+
+        { short='m' ; long="multiplicity" ; opt=Optional;
+          arg=With_arg "<int>";
+          doc="Spin multiplicity (2S+1). Default is singlet"; } ;
+
+        { short='c' ; long="charge" ; opt=Optional;
+          arg=With_arg "<int>";
+          doc="Total charge of the molecule. Default is 0"; } ;
+        
+        { short='o' ; long="output" ; opt=Mandatory;
+          arg=With_arg "<string>";
+          doc="Name of the file containing the localized MOs"; } ;
+
+        { short='i' ; long="import" ; opt=Optional;
+          arg=With_arg "<string>";
+          doc="Name of the EZFIO directory containing MOs"; } ;
+
+      ]
+  end;
+
+
+  let basis_file  = Util.of_some @@ Command_line.get "basis" in
+
+  let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
+
+  let ezfio_file = Command_line.get "import" in
+
+  let charge = 
+    match Command_line.get "charge" with
+    | Some x -> int_of_string x
+    | None   -> 0
+  in
+
+  let multiplicity =
+    match Command_line.get "multiplicity" with
+    | Some x -> int_of_string x
+    | None -> 1
+  in       
+
+  let output_filename = Util.of_some @@ Command_line.get "output" in
+
+
+  (* II : Hartree-Fock *)
+
+  (* 1. Def pour HF *)
+
+  let simulation = 
+    Simulation.of_filenames
+      ~charge ~multiplicity ~nuclei:nuclei_file basis_file 
+  in 
+
+  (* 2. Calcul de Hartree-Fock*)
+  let mo_basis = 
+     match ezfio_file with
+     | Some ezfio_file ->
+           let mo_coef =  read_qp_mo  ezfio_file in
+           let mo_type = MOBasis.Localized "Boys" in
+           let elec = Simulation.electrons simulation in
+           let n_mo = Mat.dim2 mo_coef in
+           let mo_occupation = Vec.init n_mo (fun i ->
+                   if i <= Electrons.n_beta elec then 2.
+                   else if i <= Electrons.n_alfa elec then 1.
+                   else 0.) in
+           MOBasis.make ~simulation ~mo_type ~mo_occupation ~mo_coef ()
+     | None -> HartreeFock.make simulation 
+               |> MOBasis.of_hartree_fock
+  in
+
+  let local_mos = Localisation.localize mo_basis in
+
+  let oc = open_out output_filename in
+  Printf.fprintf oc "[\n";
+  Mat.as_vec local_mos
+  |> Vec.iter (fun x -> Printf.fprintf oc "%20.15e,\n" x);
+  Printf.fprintf oc "]\n";
+  close_out oc