let png_image = print_endline ;;

(* --------- *)
    




(* --------- *)

open Lacaml.D



let charge = 0

let xyz_string = 
"3
Water
O       0.                     0.   0.
H      -0.756950272703377558   0.  -0.585882234512562827
H       0.756950272703377558   0.  -0.585882234512562827
"

(* Fonction création chaîne linéaire de n H *)

let xyz d n = 
    let accu = ""
    in
    let rec toto accu d n =
        let accu = 
            if n=0 
            then accu ^ ""
            else
                match n with 
                | 1 -> "H  0. 0. 0.\n" ^ accu
                | x -> toto ("H" ^ "  " ^ string_of_float( d *. float_of_int(n-1)) ^ "  0. 0.\n" ^ accu ) d (n-1)
            in
            accu
        in string_of_int(n) ^ "\nH" ^ string_of_int(n) ^ "\n" ^ toto accu d n;;
    
let xyz_string = xyz 1.8 6;;

let charge=1;;

let xyz_string = 
"28
Cyanine C11
C           0.000000    0.000000    0.553396
H           0.000000    0.000000    1.639052
C           0.000000    1.212451   -0.115894
C           0.000000   -1.212451   -0.115894
H           0.000000    1.184500   -1.203671
H           0.000000   -1.184500   -1.203671
C           0.000000    2.463852    0.491562
C           0.000000   -2.463852    0.491562
H           0.000000    2.518814    1.575709
H           0.000000   -2.518814    1.575709
C           0.000000    3.632731   -0.239322
C           0.000000   -3.632731   -0.239322
H           0.000000    3.543995   -1.323893
H           0.000000   -3.543995   -1.323893
C           0.000000    4.922617    0.294909
C           0.000000   -4.922617    0.294909
H           0.000000    5.053656    1.372148
H           0.000000   -5.053656    1.372148
C           0.000000    6.023589   -0.522235
C           0.000000   -6.023589   -0.522235
H           0.000000    5.885105   -1.598471
H           0.000000   -5.885105   -1.598471
N           0.000000    7.289195   -0.119132
N           0.000000   -7.289195   -0.119132
H           0.000000    7.528514    0.857192
H           0.000000   -7.528514    0.857192
H           0.000000    8.044630   -0.778923
H           0.000000   -8.044630   -0.778923
"

let xyz_string = 
"8
C1
C   0.000000    0.000000    0.424165
H   0.000000    0.000000    1.508704
N   0.000000    1.158087   -0.174215
N   0.000000   -1.158087   -0.174215
H   0.000000    1.258360   -1.178487
H   0.000000    2.005112    0.371143
H   0.000000   -2.005112    0.371143
H   0.000000   -1.258360   -1.178487
"


let basis_string =  
"
HYDROGEN
S   4
1         1.301000E+01           1.968500E-02
2         1.962000E+00           1.379770E-01
3         4.446000E-01           4.781480E-01
4         1.220000E-01           5.012400E-01
S   1
1         1.220000E-01           1.000000E+00
P   1
1         7.270000E-01           1.0000000

CARBON
S   9
1         6.665000E+03           6.920000E-04
2         1.000000E+03           5.329000E-03
3         2.280000E+02           2.707700E-02
4         6.471000E+01           1.017180E-01
5         2.106000E+01           2.747400E-01
6         7.495000E+00           4.485640E-01
7         2.797000E+00           2.850740E-01
8         5.215000E-01           1.520400E-02
9         1.596000E-01          -3.191000E-03
S   9
1         6.665000E+03          -1.460000E-04
2         1.000000E+03          -1.154000E-03
3         2.280000E+02          -5.725000E-03
4         6.471000E+01          -2.331200E-02
5         2.106000E+01          -6.395500E-02
6         7.495000E+00          -1.499810E-01
7         2.797000E+00          -1.272620E-01
8         5.215000E-01           5.445290E-01
9         1.596000E-01           5.804960E-01
S   1
1         1.596000E-01           1.000000E+00
P   4
1         9.439000E+00           3.810900E-02
2         2.002000E+00           2.094800E-01
3         5.456000E-01           5.085570E-01
4         1.517000E-01           4.688420E-01
P   1
1         1.517000E-01           1.000000E+00
D   1
1         5.500000E-01           1.0000000

NITROGEN
S   9
1         9.046000E+03           7.000000E-04
2         1.357000E+03           5.389000E-03
3         3.093000E+02           2.740600E-02
4         8.773000E+01           1.032070E-01
5         2.856000E+01           2.787230E-01
6         1.021000E+01           4.485400E-01
7         3.838000E+00           2.782380E-01
8         7.466000E-01           1.544000E-02
9         2.248000E-01          -2.864000E-03
S   9
1         9.046000E+03          -1.530000E-04
2         1.357000E+03          -1.208000E-03
3         3.093000E+02          -5.992000E-03
4         8.773000E+01          -2.454400E-02
5         2.856000E+01          -6.745900E-02
6         1.021000E+01          -1.580780E-01
7         3.838000E+00          -1.218310E-01
8         7.466000E-01           5.490030E-01
9         2.248000E-01           5.788150E-01
S   1
1         2.248000E-01           1.000000E+00
P   4
1         1.355000E+01           3.991900E-02
2         2.917000E+00           2.171690E-01
3         7.973000E-01           5.103190E-01
4         2.185000E-01           4.622140E-01
P   1
1         2.185000E-01           1.000000E+00
D   1
1         8.170000E-01           1.0000000

"

let nuclei =
  Nuclei.of_xyz_string xyz_string
  
let basis = 
  Basis.of_nuclei_and_basis_string nuclei basis_string
  
let simulation = 
  Simulation.make ~charge ~multiplicity:1 ~nuclei basis
  
let ao_basis = 
  Simulation.ao_basis simulation
  
let nocc =
    let elec = Simulation.electrons simulation in
    Electrons.n_alfa elec




let hf = HartreeFock.make simulation ;;

let ppf = Printing.ppf_dev_null ;;
Format.fprintf ppf "@[%a@]@." HartreeFock.pp hf ;;

let mo_basis = MOBasis.of_hartree_fock hf

let mo_coef = MOBasis.mo_coef mo_basis


(*
let m_X =                                                                                   
    AOBasis.ortho ao_basis
*)

(*
let c_of_h  m_H = 
    (* On exprime H dans la base orthonormale *)
    let m_Hmo =
        Util.xt_o_x m_H m_X     (* H_mo = X^t H X *)
    in
    
    (* On diagonalise cet Hamiltonien *)
    let m_C', _ =
        Util.diagonalize_symm  m_Hmo
    in
    
    (* On re-exprime les MOs dans la base des AOs (non-orthonormales) *) 
    gemm  m_X  m_C'       (* C = X.C' *)
  
  
let m_C =   
    match Guess.make ~nocc ~guess:`Hcore ao_basis with
    | Hcore m_H -> c_of_h m_H
    | _ -> assert false
*)

(*
let m_P = 
  (*  P = 2 C.C^t *)
  gemm ~alpha:2. ~transb:`T ~k:nocc  m_C m_C
  *)

(*
let m_Hc, m_J, m_K =
    let f =
      Fock.make_rhf ~density:m_P  ao_basis
    in
    Fock.(core f, coulomb f, exchange f)
    *)

(*
let m_F = 
    Mat.add m_Hc (Mat.sub m_J  m_K)
 *) 

(*let m_C =
  let m_C', _ = 
      Util.xt_o_x m_F m_X
      |> Util.diagonalize_symm
  in
  gemm  m_X  m_C'
*)

(*
let energy =
   (Simulation.nuclear_repulsion simulation) +. 0.5 *.
   Mat.gemm_trace  m_P  (Mat.add  m_Hc  m_F)
*)

(*
let rec iteration m_C n =
    let m_P = 
      (*  P = 2 C.C^t *)
      gemm ~alpha:2. ~transb:`T ~k:nocc  m_C m_C
    in

    let m_Hc, m_J, m_K =
        let f =
          Fock.make_rhf ~density:m_P  ao_basis
        in
        Fock.(core f, coulomb f, exchange f)
    in

    let m_F = 
        Mat.add m_Hc (Mat.sub m_J  m_K)
    in

    let m_C =
      let m_C', _ = 
          Util.xt_o_x m_F m_X
          |> Util.diagonalize_symm
      in
      gemm  m_X  m_C'
    in

    let energy =
       (Simulation.nuclear_repulsion simulation) +. 0.5 *.
       Mat.gemm_trace  m_P  (Mat.add  m_Hc  m_F)
    in
    Printf.printf "%f\n%!" energy;
    if n > 0 then
        iteration m_C (n-1)

let () = iteration m_C 20
*)

(* Construction de la matrice de rotation R de taille n par n *)


(* Définitions de base nécessaire pour la suite *)
let ee_ints = AOBasis.ee_ints ao_basis;;
let m_C = MOBasis.mo_coef mo_basis;;
let n_ao = Mat.dim1 m_C ;;
let n_mo = Mat.dim2 m_C ;;
let multipoles = 
      AOBasis.multipole ao_basis;;
      
let sum a = 
  Array.fold_left (fun accu x -> accu +. x) 0. a
    
let pi = 3.14159265358979323846264338;;




(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)
let f_alpha m_C =

    let n_mo = Mat.dim2 m_C in
    (*let t0 = Sys.time () 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);
   
    (*let t1 = Sys.time () in
    Printf.printf "t = %f s\n%!" (t1 -. t0);*)
    
    (Mat.init_cols n_mo n_mo ( fun i j ->
    if i= j  then 0.
    else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2. ))))
    ),Vec.sum v_d);;

(*********************)
(*
f_alpha m_C;;
let m_alpha , s_D = f_alpha m_C;;
*)

(* Fonction calcul alpha Boys *)
let f_alpha_boys m_C = 
    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 0.25 *. acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) ))),
        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.)));;

(****************************)
(*
f_alpha_boys m_C;;
*)

(* Test méthode de calcul de alpha et de D ensemble *)
type alphad = {
m_alpha : Mat.t;
d : float;
}

let m_alpha_d methode m_C = 
    let alpha methode =
        match methode with 
            | "Boys"
            | "boys" -> let alpha_boys , d_boys  = f_alpha_boys m_C in
                        {m_alpha = alpha_boys; d = d_boys}
            | "ER"
            | "er" -> let alpha_er , d_er = f_alpha m_C in
                      {m_alpha = alpha_er; d = d_er}
            | _ -> invalid_arg "Unknown method, please enter Boys or ER"

in 
alpha methode;;

(*************************)
(*
m_alpha_d "ER" ;;
m_alpha_d "Boys" ;;
let methode = "ER";;
let alphad = m_alpha_d methode m_C;;
let m_alpha = alphad.m_alpha;;
let d = alphad.d;;
*)

(* Test norme de alpha *)

let norme m = 
    let vec_m = Mat.as_vec m
    in
    let vec2 = Vec.sqr vec_m
in sqrt(Vec.sum vec2);;
 
(************************) 
 (*
norme_alpha m_alpha;;
*)

type alphaij = {
    alpha_max : float;
    indice_ii : int;
    indice_jj : int;};;

(* Détermination 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 *)
let rec new_m_alpha m_alpha m_C n_rec_alpha=

    let n_mo = Mat.dim2 m_C
    in
    let alpha_m =
        
    (*Printf.printf "%i\n%!" n_rec_alpha;*)
        
    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 if m_alpha.{i,j} < 0. 
                    then -. m_alpha.{i,j}
                    else m_alpha.{i,j} )
    in 
        
    (*Util.debug_matrix "alpha_m" alpha_m;*)
        
    (* Détermination de l'emplacement du alpha max *)
    let max_element3 alpha_m = 
        Mat.as_vec alpha_m
        |> iamax
    in
     
    (* indice i et j du alpha max *)
    let indice_ii, indice_jj = 
        let max = max_element3 alpha_m 
        in
        (max - 1) mod n_mo +1, (max - 1) / n_mo +1
    in
    
    (* Valeur du alpha max*)
    let alpha alpha_m  = 
        let i = indice_ii 
        in
        let j = indice_jj 
        in
        
        (*Printf.printf "%i %i\n%!" i j;*)
        
        alpha_m.{i,j}
            
    in
    let alpha_max = alpha alpha_m 
    in
    
    (*Printf.printf "%f\n%!" alpha_max;*)
    
    if alpha_max < pi /. 2.
        then {alpha_max; indice_ii; indice_jj}
        else new_m_alpha alpha_m m_C (n_rec_alpha-1);;

(*************************)
(*
let m_alpha,d = f_alpha m_C
let alphaij = new_m_alpha m_alpha m_C 3;;
alphaij.alpha_max;;
*)

(* Fonction de pattern matching pour localiser ou délocaliser *)
let alpha_v loc_deloc alphaij = 
    let alpha_loc = alphaij.alpha_max
    in
    let alpha_deloc = alphaij.alpha_max +. (pi /. 4.)
    in
    let choice loc_deloc = 
    match loc_deloc with
        |"loc" ->  alpha_loc
        |"deloc" ->  alpha_deloc
        | _ -> invalid_arg "Unknown method, please enter loc or deloc" 
in choice loc_deloc ;;

(* Matrice de rotation 2 par 2 *)
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 )
;;
(*************************)
(*
let alpha =  alphaij.alpha_max;;  (* Fonction -> constante *)        
f_R alpha;;
*)
(*
alpha_v "deloc" alphaij;;
let alpha = (alpha_v "loc"  alphaij)  ;;
f_R alpha ;;
*)

(*Uniquement pour pouvoir tester les fonctions après cette cellules*)
(*
(* Indice i et j du alpha max après calcul *)
let indice_i = alphaij.indice_ii;; 
let indice_j = alphaij.indice_jj;;   
let m_R = f_R alpha;;  
*)

(* 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 *)
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} );;

(*************************)
(*
let m_Ksi = f_Ksi indice_i indice_j m_C;; 
*)

(* 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*)
let f_Ksi_tilde m_R m_Ksi m_C = gemm m_Ksi m_R;;

(*************************)
(*
let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi;; 
*)

(* Fonction pour la création de matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)
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.)
(*************************)
(*
let m_Psi = f_Psi m_Ksi indice_i indice_j;; 
let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j;; 
*)

(* 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 *)
let f_interm m_C m_Psi = Mat.sub m_C m_Psi;;

(*************************)
(*
let m_interm = f_interm m_C m_Psi;; 
let new_m_C m_C= Mat.add m_Psi_tilde m_interm;;
let m_new_m_C = new_m_C m_C;;
*)

(* Test localisation matrice rectangulaire et partielle *)
(*let toto = [4];;
let occ_m_C m_C toto= Mat.init_cols 4 3 (fun i j ->
    if not (List.mem j toto) 
        then m_C.{i,j}
        else 0.);;
        
let occ = occ_m_C m_C toto;;
*)

(* Localisation de Edminstion ou de Boys *)

(* Calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)
let rec localisation m_C methode loc_deloc epsilon n prev_critere_D cc=

    (*Printf.printf "%i\n%!" n;*)

    (*Util.debug_matrix "m_C" m_C;*)

if n == 0 
    then m_C
    else
    
    (* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)
    let new_m_C m_C methode loc_deloc =
       
        (* Fonction de pattern matching en fonction de la méthode *)
        let alphad = m_alpha_d methode m_C 
        in
        
        (* D critère à maximiser *)
        let critere_D = alphad.d  
        in
  
        (*Printf.printf "%f\n%!" critere_D;*)
        
        (* Matrice des alphas *)
        let m_alpha = alphad.m_alpha
        in
        let norme_alpha = norme m_alpha
        in
        
        (*Util.debug_matrix "m_alpha" m_alpha;*)

        (* alphaij contient le alpha max ainsi que ses indices i et j *)
        let n_rec_alpha = 10 (* Nombre ditération max pour réduire les valeurs de alpha *)
        in
        let alphaij = new_m_alpha m_alpha m_C n_rec_alpha 
        in

        (* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)
        let alpha = (alpha_v loc_deloc alphaij) *. epsilon 
        in
        Printf.printf "%i %f %f %f\n%!" n critere_D alpha norme_alpha;
    

        (*Printf.printf "%f\n%!" alpha;*)
        
        (* Indice i et j du alpha max après calcul *)
        let indice_i = alphaij.indice_ii 
        in
        let indice_j = alphaij.indice_jj 
        in

        (*Printf.printf "%i %i\n%!" indice_i indice_j;*)
        
        (* Matrice de rotation *)
        let m_R = f_R alpha 
        in

        (*Util.debug_matrix "m_R" m_R;*)

        (* Matrice qui va subir la rotation *)
        let m_Ksi = f_Ksi indice_i indice_j m_C 
        in

        (*Util.debug_matrix "m_Ksi" m_Ksi;*)

        (* Matrice ayant subit la rotation *)
        let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi m_C 
        in

        (*Util.debug_matrix "m_Ksi_tilde" m_Ksi_tilde;*)

        (* Matrice pour supprimerles coef des orbitales i et j dans la matrice des coef *)
        let m_Psi = f_k m_Ksi indice_i indice_j m_C
        in

        (*Util.debug_matrix "m_Psi" m_Psi;*)

        (* Matrice pour ajouter les coef des orbitales i~ et j~ dans la matrice des coef *)
        let m_Psi_tilde = f_k m_Ksi_tilde indice_i indice_j m_C 
        in

        (*Util.debug_matrix "m_Psi_tilde" m_Psi_tilde;*)

        (* Matrice avec les coef des orbitales i et j remplacés par 0 *)
        let m_interm = f_interm m_C m_Psi 
    in
    
     (*Util.debug_matrix "m_interm" m_interm;*)
    
    (* Matrice après rotation *)
    ( Mat.add m_Psi_tilde m_interm, critere_D, norme_alpha, alpha)
    in
    let m_new_m_C , critere_D, norme_alpha, alpha = new_m_C m_C  methode loc_deloc 
    in
    let _diff = prev_critere_D -. critere_D in

    (*Util.debug_matrix "new_alpha_m" (f_alpha m_C);*)
    (*Util.debug_matrix "m_new_m_C" m_new_m_C;*)

    if alpha**2. < cc**2.
        then m_new_m_C
    else
        localisation m_new_m_C methode loc_deloc epsilon (n-1) critere_D cc;;

(* Calcul *)
(* Fonction / Matrice des coef / Méthode("Boys" ou "ER") / Localisation ou non ("loc" ou "deloc") / Pas(<=1.) 
/ Nombre d'itérations max / 0. (valeur de D pour initier la boucle) / critère de convergence sur D*)
(*
let new_m_boys = localisation m_C "boys" "loc" 1. 100  0. 10e-7;;
let new_m_er = localisation m_C "ER" "loc" 1. 100  0. 10e-7;;
*)


(*Fonction de création d'une list d'entier à partir d'un vecteur de float*)
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;;

(* Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)
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;;

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

(* Liste des OMs occupées *)
let list_occ = 
    let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))
in int_list vec_occ;;

(* Fonction de rassemblement de 2 matrices *)
let assemble m_occ m_vir =
  let occ = Mat.to_col_vecs m_occ in
  let vir = Mat.to_col_vecs m_vir in
  Array.concat [ occ ; vir ]
  |> Mat.of_col_vecs

(**********************************)
(*
m_C;;

let list_om = List.init nocc (fun i -> i+1)

let m_occ , m_vir = split_mat m_C list_om;;

assemble m_occ m_vir;;
*)


let m_occ , m_vir = split_mat m_C list_occ;;

(*
let new_m_boys = localisation m_C "boys" "loc" 1. 1000  0. 10e-7;;
let new_m_er = localisation m_C "ER" "loc" 1. 1000  0. 10e-7;;
*)

Printf.printf "Boys\n"
(*
let loc_m_occ_boys = localisation m_occ "boys" "loc" 1. 5000  0. 1e-3;;
let loc_m_vir_boys = localisation m_vir "boys" "loc" 1. 5000  0. 1e-3;;
let m_assemble_boys = assemble loc_m_occ_boys loc_m_vir_boys;;
*)

(*
let loc_m_occ_boys = localisation m_occ "boys" "loc" 1. 5000  0. 1e-3;;
let m_assemble_boys = assemble loc_m_occ_boys m_vir;;
*)

let loc_m_vir_boys = localisation m_vir "boys" "loc" 1. 50000  0. 1e-3;;
let m_assemble_boys = assemble m_occ loc_m_vir_boys;;

Printf.printf "[\n";;
Mat.as_vec m_assemble_boys
|> Vec.iter (fun x -> Printf.printf "%20.15e,\n" x);;
Printf.printf "]\n";;

(*
Printf.printf "ER\n"
let loc_m_occ_er = localisation m_occ "ER" "loc" 1. 100  0. 1e-3;;
let loc_m_vir_er = localisation m_vir "ER" "loc" 1. 100  0. 1e-3;;
let m_assemble_er = assemble loc_m_occ_er loc_m_vir_er;;

Mat.as_vec m_assemble_er
|> Vec.iter (fun x -> Printf.printf "%20.15e\n" x);;
*)

(*
let mo_base1 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized "Boys") 
~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:new_m_boys ();;

let mo_base2 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized "ER") 
~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:new_m_er ();;

let mo_base3 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized "Boys")
~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:m_assemble_boys ();;

let mo_base4 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized "ER")
~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:m_assemble_er ();;
*)

(*

(*let mo_basis = MOBasis.of_hartree_fock hf*)
let mo_basis = mo_base4

let ci =
  DeterminantSpace.fci_of_mo_basis mo_basis ~frozen_core:false
  |> CI.make
  
let ci_coef, ci_energy = Lazy.force ci.eigensystem 
let m_ci_coef = Mat.as_vec ci_coef;;
Vec.iteri (fun i x -> Printf.printf "%d %e\n" i x) m_ci_coef;;
*)