diff --git a/MOBasis/Localisation.ml b/MOBasis/Localisation.ml
index 9ddc004..a40723a 100644
--- a/MOBasis/Localisation.ml
+++ b/MOBasis/Localisation.ml
@@ -67,7 +67,7 @@
     (* 5. Fonction d'affichage *)
 
 let methode = "Boys_er";; (* Method for th orbitals localization *)
-let iteration = 10000;;   (* Maximum number of iteration for each localization *)
+let iteration = 1000000;;   (* 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. *)
@@ -79,18 +79,15 @@ let x_end = 1;; (* Choice of the x_end for the pi/sigma/end separation, size of
 (* 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 *)
+type type_separation = 
+  | Sigma_pi
+  | Occ_vir
+  | Custom of int list list
 
-(* 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 *)
+let string_of_separation = function
+  | Sigma_pi -> "Sigma/pi"
+  | Occ_vir  -> "Occ/Vir"
+  | Custom _ -> "Custom"
 
 open Lacaml.D
 
@@ -107,10 +104,10 @@ type alphad = {
   d : float;
 }
 
-let localize mo_basis =
+let localize mo_basis separation =
 
   let simulation = MOBasis.simulation mo_basis in
-  let nocc =
+  let n_occ =
       let elec = Simulation.electrons simulation
           in
           Electrons.n_alfa elec
@@ -275,52 +272,54 @@ let localize mo_basis =
       let multipoles = AOBasis.multipole ao_basis in
       let n_mo = Mat.dim2 m_C
       in
-      let phi_x_phi =
+      let x =
           Multipole.matrix_x multipoles
           |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
       in
-      let phi_y_phi =
+      let y =
           Multipole.matrix_y multipoles
           |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
       in
-      let phi_z_phi =
+      let z =
           Multipole.matrix_z multipoles
           |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
       in
-      let phi_x2_phi =
-          Multipole.matrix_x2 multipoles
+      let r2 =
+          Multipole.matrix_r2 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))
+          Mat.init_cols n_mo n_mo (fun i j ->
+              r2.{i,i} +. r2.{j,j} -.
+              (x.{i,i} *. x.{i,i} +. x.{j,j} *. x.{j,j} +. 
+               y.{i,i} *. y.{i,i} +. y.{j,j} *. y.{j,j} +.
+               z.{i,i} *. z.{i,i} +. z.{j,j} *. z.{j,j} )
+          )
       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))
+          Mat.init_cols n_mo n_mo (fun i j ->
+              let x' = (x.{i,i} -.  x.{j,j})
+              and y' = (y.{i,i} -.  y.{j,j})
+              and z' = (z.{i,i} -.  z.{j,j})
+              in
+              x'*.x' +.  y'*.y' +.  z'*.z' -.
+              4. *. ( x.{i,j} *. x.{i,j} +.
+                      y.{i,j} *. y.{i,j} +.
+                      z.{i,j} *. z.{i,j} )
+              )
       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))
+          Mat.init_cols n_mo n_mo (fun i j ->
+              4. *. ( x.{i,j} *. (x.{i,i} -. x.{j,j}) +.
+                      y.{i,j} *. (y.{i,i} -. y.{j,j}) +.
+                      z.{i,j} *. (z.{i,i} -. z.{j,j}) )
+              )
       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})
+          if i <> j then
+              0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j}
+          else 0.
+            )
       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})))
@@ -393,8 +392,8 @@ let localize mo_basis =
             then 0.
         else
           let x = atan2 m_b12.{i,j} m_a12.{i,j} in
-          if x >= 0.  then 0.25 *. (pi -. x)
-          else -. 0.25 *. (pi +. x)
+          if x >= 0. then 0.25 *. (pi -. x)
+          else         -. 0.25 *. (pi +. x)
         ),
         Vec.sum @@ Vec.init n_mo ( fun i ->
            x.{i,i} *. x.{i,i} +.
@@ -466,7 +465,7 @@ let localize mo_basis =
                 else +. 0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j})
         in
         let m_critere_B = Vec.init n_mo (fun i ->
-            r2.{i,i} -. x.{i,i}**2. -. y.{i,i}**2. -. z.{i,i}**2.)
+            r2.{i,i} -. x.{i,i}*.x.{i,i} -. y.{i,i}*.y.{i,i} -. z.{i,i}*.z.{i,i})
 
         in m_theta, Vec.sum( m_critere_B)
    in
@@ -512,11 +511,7 @@ let localize mo_basis =
   (* 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)
+      sqrt (Mat.syrk_trace m)
   in
 
   (* 7. Calculation of the max angle[-pi/.4 , pi/.4]  and its indexes i and j in the matrix *)
@@ -577,13 +572,9 @@ let localize mo_basis =
   (* 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 )
+      [| [|  cos alpha ; -. sin alpha |] ; 
+         [|  sin alpha ;    cos alpha |] |]
+      |> Mat.of_array 
   in
 
   (* 9. Rotation matrix nxn *)
@@ -902,28 +893,6 @@ let localize mo_basis =
   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 lst =
-
-      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 lst
-      in
-      Mat.of_col_vecs_list vec_list
-  in
-
-  (* List of the occupied MOs *)
-  let list_occ =
-    Util.list_range 1 nocc 
-  in
-
-  (* List of the core MOs *)
-  let list_core =
-    Util.list_range 1 n_core
-  in
-
   (* Function to create a list of pi MOs *)
   (* Coefficient matrix n_ao x n_mo -> list of pi orbitals *)
   (* matrix -> integer list *)
@@ -933,17 +902,17 @@ let localize mo_basis =
           if j=1 && abs_float m_C.{i,j} < 1.e-8
               then 0.
           else if j=1
-              then m_C.{i,j}
+              then abs_float m_C.{i,j}
           else if abs_float mat.{i,j-1} < 1.e-8
               then 0.
-              else mat.{i,j-1})
+              else abs_float mat.{i,j-1})
       in
       let vec = Mat.to_col_vecs m_pi in
       (* TODO gemv *)
       let vec_dot = Vec.init (Mat.dim2 mat) (fun i -> dot vec.(0) vec.(i)) in
       let list_pi =
         Util.list_range 1 (Mat.dim2 mat) 
-        |> List.map (fun i -> if vec_dot.{i} = 0. then i else 0)
+        |> List.map (fun i -> if abs_float vec_dot.{i} = 0. then i else 0)
       in
       let rec remove x = function
         | [] -> []
@@ -954,22 +923,6 @@ let localize mo_basis =
       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
-      Util.list_range 1 (n-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
-      Util.list_range 1 (n/2)
-  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*)
@@ -985,23 +938,22 @@ let localize mo_basis =
 (* 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)
+  let text = "Method : " ^ methode ^ "\nSeparation type : "
+    ^ (string_of_separation separation)   
+    ^ "\nMax number of iterations : " ^ string_of_int(iteration) ^ "\ncc : "
+    ^ string_of_float(cc)^"\nepsilon : "^string_of_float(epsilon)^"\ncharge : "
+    ^ string_of_int(charge)^"\nRotation pre angle : "^string_of_float(pre_angle)
   in
-  let caption = "n Criterion max_angle norm_angle average_angle"
+  let caption = "Criterion max_angle norm_angle average_angle"
   in
   Printf.printf "%s\n" "";
-  Printf.printf "%s\n" "Orbitals localization ";
+  Printf.printf "%s\n" "Orbital 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 *)
@@ -1020,215 +972,68 @@ let localize mo_basis =
           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
+  (* Localization function with n separations expressed as a list of lists of
+     orbital indices *)
+  let f_loc_list mat list_of_lists =
+    let vec_list = Mat.to_col_vecs mat in
+    list_of_lists
+    |> List.map (fun lst ->
+           lst
+           |> List.map (fun col -> vec_list.(col-1)) 
+           |> Mat.of_col_vecs_list
+           |> f_localise
+           |> Mat.to_col_vecs_list
+         )
+    |> List.concat
+    |> Mat.of_col_vecs_list
   in
 
-  (* Localization function with 2 MOs separations *)
-  let f_loc_assemble_123 mat list_12 list_2ab =
+  let n_mo = Mat.dim2 m_C in
+  let list_full      = Util.list_range 1 n_mo in
+  let list_core      = Util.list_range 1 n_core in
+  let list_pi        = list_pi m_C in
+  let list_sigma     = List.filter (fun i -> not (List.mem i list_pi)) list_full in
+  let list_pi_occ    = List.filter (fun i -> i > n_core && i <= n_occ) list_pi in
+  let list_sigma_occ = List.filter (fun i -> i > n_core && i <= n_occ) list_sigma in
+  let list_pi_vir    = List.filter (fun i -> i > n_occ) list_pi in
+  let list_sigma_vir = List.filter (fun i -> i > n_occ) list_sigma in
+  let list_occ       = Util.list_range (n_core+1) n_occ in
+  let list_vir       = Util.list_range (n_occ+1) n_mo in
 
-      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)";
+  Printf.printf "Orbital lists:";
+  Printf.printf "\nCore: [";
+  List.iter (fun i -> Printf.printf "%d," i) list_core;
+  Printf.printf "\b]\nOccupied: [";
+  List.iter (fun i -> Printf.printf "%d," i) list_occ;
+  Printf.printf "\b]\nVirtual: [";
+  List.iter (fun i -> Printf.printf "%d," i) list_vir;
+  Printf.printf "\b]\nSigma occ: [";
+  List.iter (fun i -> Printf.printf "%d," i) list_sigma_occ;
+  Printf.printf "\b]\nPi occ: [";
+  List.iter (fun i -> Printf.printf "%d," i) list_pi_occ;
+  Printf.printf "\b]\nSigma vir: [";
+  List.iter (fun i -> Printf.printf "%d," i) list_sigma_vir;
+  Printf.printf "\b]\nPi vir: [";
+  List.iter (fun i -> Printf.printf "%d," i) list_pi_vir;
+  Printf.printf "\b]\n\n";
 
-      let mat_1, mat_2 = split_mat mat (list_12 mat) in
+  (* 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
 
-      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
+  let m_assemble_loc = 
+    let list_of_lists = 
+      match separation with
+        | Occ_vir        -> [ list_core ; list_occ ; list_vir ]
+        | Sigma_pi       -> [ list_core ; list_sigma_occ ; list_pi_occ ; list_pi_vir ; list_sigma_vir ]
+        | Custom l       -> l
+    in
+    f_loc_list m_C list_of_lists
   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
+  let m_assemble_occ, m_assemble_vir = split_mat m_assemble_loc list_occ in
+  let m_assemble_core, m_assemble_occu = split_mat m_assemble_occ list_core in
 
   (* VI. Analysis *)
 
@@ -1337,6 +1142,7 @@ let localize mo_basis =
       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 *)
@@ -1345,19 +1151,20 @@ let localize mo_basis =
   (* 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 %i\n" "Number of pi orbitals : " (List.length list_pi );
+  Printf.printf "%s %i\n" "Number of occupied pi orbitals : " (List.length list_pi_occ );
+  Printf.printf "%s %i\n" "Number of virtual pi orbitals : " (List.length list_pi_vir );
+  Printf.printf "%s %i\n" "Number of sigma orbitals : " (List.length list_sigma);
+  Printf.printf "%s %i\n" "Number of occupied sigma orbitals : " (List.length list_sigma_occ);
+  Printf.printf "%s %i\n" "Number of virtual sigma orbitals : " (List.length list_sigma_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));
+  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);
@@ -1368,7 +1175,7 @@ let localize mo_basis =
   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 %s\n" "Core orbitals after localization : " (analyse m_assemble_core);
   Printf.printf "%s\n" "";
 
   Printf.printf "%s\n" "";
diff --git a/run_Localisation.ml b/run_Localisation.ml
index 32eb5ae..e6679f7 100644
--- a/run_Localisation.ml
+++ b/run_Localisation.ml
@@ -63,6 +63,10 @@ let () =
           arg=With_arg "<string>";
           doc="Name of the EZFIO directory containing MOs"; } ;
 
+        { short='s' ; long="separation" ; opt=Optional;
+          arg=With_arg "<string>";
+          doc="Type of separation: [ov: occ/vir | sp: sigma/pi | custom] "; } ;
+
       ]
   end;
 
@@ -87,6 +91,29 @@ let () =
 
   let output_filename = Util.of_some @@ Command_line.get "output" in
 
+  let separation =
+    match Command_line.get "separation" with
+      | None
+      | Some "ov" -> Localisation.Occ_vir
+      | Some "sp" -> Localisation.Sigma_pi
+      | Some "custom" ->
+        let list_of_lists =
+          print_endline "Enter ranges, one per line:";
+          let rec loop accu =
+            try
+              let r =
+                input_line stdin
+                |> Range.of_string
+                |> Range.to_int_list
+              in
+              loop (r::accu)
+            with End_of_file -> List.rev accu
+          in loop []
+        in
+        Localisation.Custom list_of_lists
+      | _ -> invalid_arg "separation"
+  in
+    
 
   (* II : Hartree-Fock *)
 
@@ -114,7 +141,7 @@ let () =
                |> MOBasis.of_hartree_fock
   in
 
-  let local_mos = Localisation.localize mo_basis in
+  let local_mos = Localisation.localize mo_basis separation in
 
   let oc = open_out output_filename in
   Printf.fprintf oc "[\n";