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) = = \sum_q 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) = = \sum_r 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;; +*) + + + + + +