(* 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 = 1_000_000;;   (* 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 *)
type type_separation = 
  | Sigma_pi
  | Occ_vir
  | Custom of int list list

let string_of_separation = function
  | Sigma_pi -> "Sigma/pi"
  | Occ_vir  -> "Occ/Vir"
  | Custom _ -> "Custom"

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 separation =

  let simulation = MOBasis.simulation mo_basis in
  let n_occ =
      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
    *)
    0
  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 x =
          Multipole.matrix_x multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
      in
      let y =
          Multipole.matrix_y multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
      in
      let z =
          Multipole.matrix_z multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
      in
      let m_b12=
          Mat.init_cols n_mo n_mo ( fun i j ->
            (x.{i,i} -. x.{j,j}) *. x.{i,j} +.
            (y.{i,i} -. y.{j,j}) *. y.{i,j} +.
            (z.{i,i} -. z.{j,j}) *. z.{i,j})
      in
      let m_a12 =
          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.{i,j} *. x.{i,j} +. 
              y.{i,j} *. y.{i,j} +. 
              z.{i,j} *. z.{i,j} -. 
              0.25 *. (x'*.x' +. y'*.y' +. z'*.z')
            )
      in
      (Mat.init_cols n_mo n_mo ( fun i j ->
          if i=j
              then 0.
          else
            let w = atan2 m_b12.{i,j} m_a12.{i,j} in
            if  w >= 0.
            then    0.25 *. (pi  -. w)
            else -. 0.25 *. (pi +. w) 
        ) , 0.5 *. (Vec.sum @@ Vec.init n_mo ( fun i ->
               x.{i,i} *. x.{i,i} +.
               y.{i,i} *. y.{i,i} +. 
               z.{i,i} *. z.{i,i} ) )
      )
  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 x =
          Multipole.matrix_x multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
      in
      let y =
          Multipole.matrix_y multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
      in
      let z =
          Multipole.matrix_z multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
      in
      let r2 =
          Multipole.matrix_r2 multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
      in
      let m_b_0 =
          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 =
          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 =
          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.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})))

      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 x =
          Multipole.matrix_x multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
    in
    let y =
          Multipole.matrix_y multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
    in
    let z =
          Multipole.matrix_z multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
    in

    let m_b12 = 
      Mat.init_cols n_mo 2 (fun i j ->
        if j = 1 then
          (x.{i,i} -. x.{j1,j1}) *. x.{i,j1} +. 
          (y.{i,i} -. y.{j1,j1}) *. y.{i,j1} +. 
          (z.{i,i} -. z.{j1,j1}) *. z.{i,j1}
        else
          (x.{i,i} -. x.{j2,j2}) *. x.{i,j2} +. 
          (y.{i,i} -. y.{j2,j2}) *. y.{i,j2} +. 
          (z.{i,i} -. z.{j2,j2}) *. z.{i,j2}
        )
    in

    let m_a12 =
      Mat.init_cols n_mo 2 (fun i j ->
        if j = 1 then

          let x' = (x.{i,i} -. x.{j1,j1})
          and y' = (y.{i,i} -. y.{j1,j1})
          and z' = (z.{i,i} -. z.{j1,j1}) in
          x.{i,j1} *. x.{i,j1} +.
          y.{i,j1} *. y.{i,j1} +. 
          z.{i,j1} *. z.{i,j1} -.
          0.25 *. (x' *. x' +. y' *. y' +. z' *. z' )

        else

          let x' = (x.{i,i} -. x.{j2,j2})
          and y' = (y.{i,i} -. y.{j2,j2})
          and z' = (z.{i,i} -. z.{j2,j2}) in
          x.{i,j2} *. x.{i,j2} +.
          y.{i,j2} *. y.{i,j2} +. 
          z.{i,j2} *. z.{i,j2} -.
          0.25 *. (x' *. x' +. y' *. y' +. z' *. z' )
       )
    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
          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)
        ),
        Vec.sum @@ Vec.init n_mo ( fun i ->
           x.{i,i} *. x.{i,i} +.
           y.{i,i} *. y.{i,i} +.
           z.{i,i} *. z.{i,i} )
        )
  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 x =
              Multipole.matrix_x multipoles
              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
        in
        let y =
              Multipole.matrix_y multipoles
              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
        in
        let z =
              Multipole.matrix_z multipoles
              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
        in
        let r2 =
              Multipole.matrix_r2 multipoles
              |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
        in

        let m_beta =
            Mat.init_cols n_mo 2 (fun i j ->
              if j = 1 then
           
                let x' = (x.{i,i} -. x.{j1,j1})
                and y' = (y.{i,i} -. y.{j1,j1})
                and z' = (z.{i,i} -. z.{j1,j1}) in
                x' *. x' +. y' *. y' +. z' *. z' 
                -. 4. *. ( x.{i,j1} *. x.{i,j1} +.
                           y.{i,j1} *. y.{i,j1} +. 
                           z.{i,j1} *. z.{i,j1} )
              else
           
                let x' = (x.{i,i} -. x.{j2,j2})
                and y' = (y.{i,i} -. y.{j2,j2})
                and z' = (z.{i,i} -. z.{j2,j2}) in
                x' *. x' +. y' *. y' +. z' *. z' 
                -. 4. *. ( x.{i,j2} *. x.{i,j2} +.
                           y.{i,j2} *. y.{i,j2} +. 
                           z.{i,j2} *. z.{i,j2} )
            )
        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 x) (Mat.add (gamma y) (gamma z))
        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 ->
            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

  (* 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 =
      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 *)
  (* 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 pi_2 = 0.5 *. pi in
      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  =
      [| [|  cos alpha ; -. sin alpha |] ; 
         [|  sin alpha ;    cos alpha |] |]
      |> Mat.of_array 
  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) *)
      (* -> float *)
      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 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

      (* W cos(tau) Wt *)
      (* -> matrix *)
      let a = gemm mm @@ gemm ~transb:`T m_cos_tau mm in

      (* W tau^-1 sin(tau) Wt X *)
      (* -> matrix *)
      let b = gemm mm @@ gemm m_tau_1 @@ gemm m_sin_tau @@ gemm ~transa:`T 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 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 lst =

      let n_mo = Mat.dim2 mat in
      Util.list_range 1 n_mo
      |> List.map (fun i ->
          if List.mem i lst
              then 0
              else i)
      |> List.filter (fun x -> x > 0) 
  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 lst =

      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 lst in
      let list_2 = miss_elem mat lst 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 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 && abs_float m_C.{i,j} < 1.e-8
              then 0.
          else if j=1
              then abs_float m_C.{i,j}
          else if abs_float mat.{i,j-1} < 1.e-8
              then 0.
              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 abs_float vec_dot.{i} = 0. then i else 0)
      in
      let rec remove x = function
        | [] -> []
        | var :: tail -> if var = x
                       then remove x tail
                       else var :: (remove x tail)
      in
      remove 0 list_pi
  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 ^ "\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 = "Criterion max_angle norm_angle average_angle"
  in
  Printf.printf "%s\n" "";
  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" "";

  (* 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


  (* 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

  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

  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";

  (* 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

  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

  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 *)

  (* 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 x =
          Multipole.matrix_x multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
      in
      let y =
          Multipole.matrix_y multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
      in
      let z =
          Multipole.matrix_z multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
      in
      let r2 =
          Multipole.matrix_r2 multipoles
          |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
      in
      Vec.init n_mo (fun i ->
      r2.{i,i} -. x.{i,i}*.x.{i,i} -. y.{i,i}*.y.{i,i} -. z.{i,i}*.z.{i,i}))
  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)) <= 1
          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) - 1)
  in

  (* Function to compute the standard deviation *)
  let f_stand_dev mat =
          if (Vec.dim (distrib mat)) <= 1
          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 );
  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));

  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 m_assemble_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