46 KiB
46 KiB
None
<html>
<head>
</head>
</html>
Sauvegarde des anciennes fonctions¶
In [ ]:
In [ ]:
(*
(* 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 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 =
(* 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
Printf.printf "%f\n%!" norme_alpha;
(*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 = (alphaij.alpha_max) *. epsilon
in
(*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)
in
let m_new_m_C , critere_D = new_m_C m_C methode
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 diff**2. < cc**2.
then m_new_m_C
else
localisation m_new_m_C methode epsilon (n-1) critere_D cc;;
*)
In [ ]:
(* Fonction de calcul de D Boys *)
(*
let d_boys m_C =
let phi_x_phi =
Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
(*Util.debug_matrix "phi_x_phi" phi_x_phi;*)
let phi_y_phi =
Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
(*Util.debug_matrix "phi_y_phi" phi_y_phi;*)
let phi_z_phi =
Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
(*Util.debug_matrix "phi_z_phi" phi_z_phi;*)
let v_D_boys =
let n_mo = Mat.dim2 m_C
in
Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)
in
Vec.sum v_D_boys;;
*)
(*************************)
(*
Multipole.matrix_x multipoles;;
d_boys m_C;;
*)
In [ ]:
(* 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 vec_list = Vec.to_list vec
in
let g a = int_of_float(a)
in
let vec_list_int = List.map g vec_list
*)
In [ ]:
(* 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
(* 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)
) (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. ))))
);;
(*********************)
f_alpha m_C;;
*)
In [ ]:
(*
(* 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 eps=
let n_mo = Mat.dim2 m_C
in
(* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)
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
Array.iter (fun a ->
Array.iter (fun b ->
let mca = Vec.init n_mo (fun i -> m_C.{a,i} *. m_C.{b,i}) in
Array.iter (fun e ->
Array.iter (fun f ->
let integral = ERI.get_phys ee_ints f b e a in
if abs_float integral > eps then
Array.iter ( fun i ->
let mcei = m_C.{e,i} *. m_C.{f,i} in
Array.iter ( fun j ->
let x = m_C.{e,i} *. m_C.{f,j} *. integral in
let mcaij = ( mca.{i} -. mca.{j} ) in
m_b12.{i,j} <- m_b12.{i,j} +. mcaij *. x;
m_a12.{i,j} <- m_a12.{i,j} +. m_C.{a,i} *. m_C.{b,j} *. x
-. 0.25 *. ( mcei -. m_C.{e,j} *. m_C.{f,j} ) *. mcaij *. integral
) (Util.array_range 1 n_mo)
) (Util.array_range 1 n_mo)
) (Util.array_range 1 n_ao)
) (Util.array_range 1 n_ao)
) (Util.array_range 1 n_ao)
) (Util.array_range 1 n_ao);
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.)))));;
(*********************)
f_alpha m_C 1.e-5;;
*)
In [ ]:
(*
(*Fonction général de calcul des intégrales*)
let integral_general g i j =
Array.map (fun a ->
let v =
Array.map (fun b ->
let u =
Array.map (fun e ->
let t = Array.map (fun f ->
(g a b e f i j) *. ERI.get_phys ee_ints a e b f
) (Util.array_range 1 n_ao)
in sum t
) (Util.array_range 1 n_ao)
in sum u
) (Util.array_range 1 n_ao)
in sum v
) (Util.array_range 1 n_ao)
|> sum
(* 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
(* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)
let m_b12 = Mat.init_cols n_mo n_mo (fun i j ->
integral_general (fun a b e f i j ->
( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}
) i j
)
in
(* Fonction de calcul de toutes les intégrales A_12 -> Matrice, dépend de m_C *)
let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j} *. m_C.{e,i} *. m_C.{f,j}
-. 0.25 *. (( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} )
*. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ))
) i j
)
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.)))));;
(*********************)
(*
f_alpha m_C;;
*)
*)
In [ ]:
(*
(* Calcul de D -> critère à maximiser dans ER*)
let s_D m_C =
let v_D =
let n_mo = Mat.dim2 m_C
in
let m_D = Mat.init_cols n_mo n_mo (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i}
) i j
)
in Vec.init n_mo ( fun i -> m_D.{i,i} )
in Vec.sum v_D ;;
(******************)
let m_D = Mat.init_cols n_mo n_mo (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i}
) i j
);;
let toto = s_D m_C;;
*)
Boucle localisation avec fonctions définies à l'extérieur¶
In [ ]:
(* Ancienne boucle
(* Localisation de Edminstion *)
let n_rec_alpha = 10;;
(* Calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)
let rec final_m_C m_C n =
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 =
(* D critère à maximiser *)
let critere_D = s_D m_C (* Fonction -> constante *)
in
Printf.printf "%f\n%!" critere_D;
(* Matrice des alphas *)
let m_alpha = f_alpha m_C (* Fonction -> constante *)
in
(*Util.debug_matrix "m_alpha" m_alpha;*)
(* alphaij contient le alpha max ainsi que ses indices i et j *)
let alphaij = new_m_alpha m_alpha n_rec_alpha (* Fonction -> constante *)
in
(* Valeur de alpha max après calcul *)
let alpha = alphaij.alpha_max (* Fonction -> constante *)
in
(* Indice i et j du alpha max après calcul *)
let indice_i = alphaij.indice_ii (* Fonction -> constante *)
in
let indice_j = alphaij.indice_jj (* Fonction -> constante *)
in
(*Printf.printf "%i %i\n%!" indice_i indice_j;*)
(* Matrice de rotaion *)
let m_R = f_R alpha (* Fonction -> constante *)
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 (* Fonction -> constante *)
in
(*Util.debug_matrix "m_Ksi" m_Ksi;*)
(* Matrice ayant subit la rotation *)
let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi (* Fonction -> constante *)
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_Psi m_Ksi indice_i indice_j (* Fonction -> constante *)
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_Psi_tilde m_Ksi_tilde indice_i indice_j (* Fonction -> constante *)
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 (* Fonction -> constante *)
in
(* Matrice après rotation *)
Mat.add m_Psi_tilde m_interm
in
let m_new_m_C = new_m_C m_C (* Fonction -> constante *)
in
(*Util.debug_matrix "new_alpha_m" (f_alpha m_C);*)
(*Util.debug_matrix "m_new_m_C" m_new_m_C;*)
final_m_C m_new_m_C (n-1);;
*)
In [ ]:
(*
(*Cellule de calcul*)
final_m_C m_C 10;;
*)
In [ ]:
(*Définition des fonctions teste*)
(*
(* Localisation de Edminstion *)
(* 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 sum a =
Array.fold_left (fun accu x -> accu +. x) 0. a
(******************************************************************************************************)
(*Fonction général de calcul des intégrales*)
let integral_general g i j =
Array.map (fun a ->
let v =
Array.map (fun b ->
let u =
Array.map (fun e ->
let t = Array.map (fun f ->
(g a b e f i j) *. ERI.get_phys ee_ints a e b f
) (Util.array_range 1 n_ao)
in sum t
) (Util.array_range 1 n_ao)
in sum u
) (Util.array_range 1 n_ao)
in sum v
) (Util.array_range 1 n_ao)
|> sum
;;
(******************************************************************************************************)
(* Fonction de calcul de tous les alpha -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)
let f_alpha m_C =
(* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)
let m_b12 = Mat.init_cols n_ao n_ao (fun i j ->
integral_general (fun a b e f i j ->
( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}
) i j
)
in
(* Fonction de calcul de toutes les intégrales A_12 -> Matrice, dépend de m_C *)
let m_a12 = Mat.init_cols n_ao n_ao (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j} *. m_C.{e,i} *. m_C.{f,j}
-. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} )
*. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )
) i j
)
in
Mat.init_cols n_ao n_ao ( 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.)))))
;;
f_alpha m_C;;
(******************************************************************************************************)
(* Calcul de D *)
let s_D m_C =
let v_D =
let m_D = Mat.init_cols n_ao n_ao (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i}
) i j
)
in Vec.init n_ao ( fun i -> m_D.{i,1} )
in Vec.sum v_D ;;
s_D m_C;;
(******************************************************************************************************)
(* Si alpha max > pi/2 on doit soustraire pi/2 à la matrice des alphas *)
let m_alpha = f_alpha m_C;;
type alphaij = {
alpha_max : float;
indice_ii : int;
indice_jj : int;}
let rec new_m_alpha m_alpha n_rec_alpha=
let alpha_m =
Printf.printf "%i\n%!" n_rec_alpha;
if n_rec_alpha == 0
then m_alpha
else Mat.init_cols n_ao n_ao (fun i j ->
if (m_alpha.{i,j}) > 3.14 /. 2.
then (m_alpha.{i,j} -. ( 3.14 /. 2.))
else if m_alpha.{i,j} < -. 3.14 /. 2.
then (m_alpha.{i,j} +. ( 3.14 /. 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;
let max_element3 alpha_m =
Mat.as_vec alpha_m
|> iamax
in
let indice_ii =
let max = max_element3 alpha_m
in
Printf.printf "%i\n%!" max;
(max - 1) mod n_ao +1
in
let indice_jj =
let max = max_element3 alpha_m
in
(max - 1) / n_ao +1
in
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
Printf.printf "%f\n%!" alpha_max;
if alpha_max < 3.14 /. 2.
then {alpha_max; indice_ii; indice_jj}
else new_m_alpha alpha_m (n_rec_alpha-1);;
new_m_alpha m_alpha 10 ;;
let alphaij = new_m_alpha m_alpha 10 ;;
(******************************************************************************************)
let alpha = alphaij.alpha_max
(* Matrice de rotation 2 par 2 *)
let m_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 indice_i = alphaij.indice_ii;;
let indice_j = alphaij.indice_jj;;
let m_R = m_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 m_Ksi indice_i indice_j m_C = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_C.{i,indice_i} else m_C.{i,indice_j} );;
m_Ksi indice_i indice_j m_C;;
let m_Ksi = m_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 m_Ksi_tilde m_R m_Ksi = gemm m_Ksi m_R;;
m_Ksi_tilde m_R m_Ksi ;;
(******************************************************************************************)
let m_Ksi_tilde = m_Ksi_tilde m_R m_Ksi
(* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur
les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice
des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec
celle contenant i~ et j~ *)
(* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)
let m_Psi_tilde m_Ksi_tilde indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi_tilde.{i,1}
else if j=indice_j then m_Ksi_tilde.{i,2}
else 0.);;
m_Psi_tilde m_Ksi_tilde indice_i indice_j;;
(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *)
let m_Psi m_Ksi indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi.{i,1}
else if j=indice_j then m_Ksi.{i,2}
else 0.);;
m_Psi m_Ksi indice_i indice_j;;
(******************************************************************************************)
let m_Psi = m_Psi m_Ksi indice_i indice_j;;
let m_Psi_tilde = m_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 m_interm m_C m_Psi = Mat.sub m_C m_Psi;;
let m_interm = m_interm m_C m_Psi;;
(* Construction de la nouvelle matrice d'OMs phi~ *)
let m_Phi_tilde m_Psi_tilde m_interm = Mat.add m_Psi_tilde m_interm;;
m_Phi_tilde m_Psi_tilde m_interm;;
*)
Programmes avec fonctions définies à l'intérieurs¶
In [ ]:
(* Localisation de Edminstion *)
(*Fonctionne -> programme avec les fonctions définies à l'intérieur*)
(*
(* 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 sum a =
Array.fold_left (fun accu x -> accu +. x) 0. a
(******************************************************************************************************)
(*Fonction général de calcul des intégrales*)
let integral_general g i j =
Array.map (fun a ->
let v =
Array.map (fun b ->
let u =
Array.map (fun e ->
let t = Array.map (fun f ->
(g a b e f i j) *. ERI.get_phys ee_ints a e b f
) (Util.array_range 1 n_ao)
in sum t
) (Util.array_range 1 n_ao)
in sum u
) (Util.array_range 1 n_ao)
in sum v
) (Util.array_range 1 n_ao)
|> sum
;;
type alphaij = {
alpha_max : float;
indice_ii : int;
indice_jj : int;};;
let n_rec_alpha = 10;;
(*
let alpha = 1.;;
let y_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 y_R = y_R alpha;;
let m_C = gemm m_C y_R;;
*)
(******************************************************************************************************)
let rec final_m_C m_C n =
Printf.printf "%i\n%!" n;
(*
Util.debug_matrix "m_C" m_C;
*)
if n == 0 then m_C
else
let new_m_C m_C =
(* Fonction de calcul de tous les alpha -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)
let f_alpha m_C =
(* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)
let m_b12 = Mat.init_cols n_ao n_ao (fun i j ->
integral_general (fun a b e f i j ->
( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}
) i j
)
in
(* Fonction de calcul de toutes les intégrales A_12 -> Matrice, dépend de m_C *)
let m_a12 = Mat.init_cols n_ao n_ao (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j} *. m_C.{e,i} *. m_C.{f,j}
-. 0.25 *. (( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} )
*. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ))
) i j
)
in
Mat.init_cols n_ao n_ao ( 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.)))))
in
(* Calcul de D *)
let s_D m_C =
let v_D =
let m_D = Mat.init_cols n_ao n_ao (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i}
) i j
)
in Vec.init n_ao ( fun i -> m_D.{i,1} )
in Vec.sum v_D
in
let critere_D = s_D m_C (* Fonction -> constante *)
in
Printf.printf "%f\n%!" critere_D;
let m_alpha = f_alpha m_C (* Fonction -> constante *)
in
(*
Util.debug_matrix "m_alpha" m_alpha;
*)
(* 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 n_rec_alpha=
let alpha_m =
Printf.printf "%i\n%!" n_rec_alpha;
if n_rec_alpha == 0
then m_alpha
else Mat.init_cols n_ao n_ao (fun i j ->
if (m_alpha.{i,j}) > (3.14 /. 2.)
then (m_alpha.{i,j} -. ( 3.14 /. 2.))
else if m_alpha.{i,j} < -. 3.14 /. 2.
then (m_alpha.{i,j} +. ( 3.14 /. 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 du alpha max *)
let indice_ii =
let max = max_element3 alpha_m (* Fonction -> constante *)
in
(*
Printf.printf "%i\n%!" max;
*)
(max - 1) mod n_ao +1
in
(* indice j du alpha max *)
let indice_jj =
let max = max_element3 alpha_m (* Fonction -> constante *)
in
(max - 1) / n_ao +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 (* Fonction -> constante *)
in
Printf.printf "%f\n%!" alpha_max;
if alpha_max < 3.14 /. 2.
then {alpha_max; indice_ii; indice_jj}
else new_m_alpha alpha_m (n_rec_alpha-1)
in
let alphaij = new_m_alpha m_alpha n_rec_alpha (* Fonction -> constante *)
in
(* Valeur de alpha max après calcul *)
let alpha = alphaij.alpha_max (* Fonction -> constante *)
in
(* 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 )
in
(* Indice i et j du alpha max après calcul *)
let indice_i = alphaij.indice_ii (* Fonction -> constante *)
in
let indice_j = alphaij.indice_jj (* Fonction -> constante *)
in
Printf.printf "%i %i\n%!" indice_i indice_j;
let m_R = f_R alpha (* Fonction -> constante *)
in
Util.debug_matrix "m_R" m_R;
(* 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 = 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
let m_Ksi = f_Ksi indice_i indice_j m_C (* Fonction -> constante *)
in
Util.debug_matrix "m_Ksi" m_Ksi;
(* 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 = gemm m_Ksi m_R
in
let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi (* Fonction -> constante *)
in
Util.debug_matrix "m_Ksi_tilde" m_Ksi_tilde;
(* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur
les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice
des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec
celle contenant i~ et j~ *)
(* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)
let f_Psi_tilde m_Ksi_tilde indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi_tilde.{i,1}
else if j=indice_j then m_Ksi_tilde.{i,2}
else 0.)
in
(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *)
let f_Psi m_Ksi indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi.{i,1}
else if j=indice_j then m_Ksi.{i,2}
else 0.)
in
let m_Psi = f_Psi m_Ksi indice_i indice_j (* Fonction -> constante *)
in
Util.debug_matrix "m_Psi" m_Psi;
let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j (* Fonction -> constante *)
in
Util.debug_matrix "m_Psi_tilde" m_Psi_tilde;
(* 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
in
let m_interm = f_interm m_C m_Psi (* Fonction -> constante *)
in
Mat.add m_Psi_tilde m_interm
in
let m_new_m_C = new_m_C m_C (* Fonction -> constante *)
in
Util.debug_matrix "new_alpha_m" (f_alpha m_C);
Util.debug_matrix "m_new_m_C" m_new_m_C;
final_m_C m_new_m_C (n-1);;
(*****************************)
final_m_C m_C 10;;
*)
In [ ]:
(*
(* Localisation de Edminstion *)
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 sum a =
Array.fold_left (fun accu x -> accu +. x) 0. a
;;
(*Calcul des intégrales*)
let integral_general g i j =
Array.map (fun a ->
let v =
Array.map (fun b ->
let u =
Array.map (fun e ->
let t = Array.map (fun f ->
(g a b e f i j) *. ERI.get_phys ee_ints a e b f
) (Util.array_range 1 n_ao)
in sum t
) (Util.array_range 1 n_ao)
in sum u
) (Util.array_range 1 n_ao)
in sum v
) (Util.array_range 1 n_ao)
|> sum
;;
(* Calcul de toutes les intégrales B_12 -> Matrice *)
let m_b12 = Mat.init_cols n_ao n_ao (fun i j ->
integral_general (fun a b e f i j ->
( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}
) i j
)
(* Calcul de toutes les intégrales A_12 -> Matrice *)
let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j} *. m_C.{e,i} *. m_C.{f,j}
-. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} )
*. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )
) i j
)
(* Calcul de tous les alpha -> Matrice *)
let m_alpha = 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.)))))
(*
let s_D =
let v_D =
let m_D = Mat.init_cols n_mo n_mo (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i}
) i j
)
in Vec.init n_mo ( fun i -> m_D.{i,1} )
in Vec.sum v_D
*)
(* Calcul de D *)
let m_D = Mat.init_cols n_mo n_mo (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i}
) i j
)
let v_D = Vec.init n_mo ( fun i -> m_D.{i,1} )
let s_D = Vec.sum v_D
(*
(* Calcul de D *)
let v_D = Vec.init n_ao (fun i ->
integral_general (fun a b e f i _ -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i}
)
)
let s_D = Vec.sum v_D
*)
(* Rotation des orbitales *)
(* Extraction du plus grand angle alpha*)
let ij_max_alpha = Mat.as_vec m_alpha
|> iamax;;
(* Indices du plus grand angle alpha *)
let indice_i = (ij_max_alpha - 1) mod n_ao + 1;;
let indice_j = (ij_max_alpha - 1) / n_mo + 1;;
(* Valeur du plus grand angle alpha*)
let alpha = m_alpha.{indice_i,indice_j};;
let m_Phi = m_C
(* Matrice de rotation 2 par 2 *)
let m_R = Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha
else if i>j then sin alpha
else -. sin 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 p = indice_i;;
let q = indice_j;;
let m_Ksi = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_Phi.{i,p} else m_Phi.{i,q} );;
(* 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 m_Ksi_tilde = gemm m_Ksi m_R;;
(* Réinjection*)
(* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur
les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice
des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec
celle contenant i~ et j~ *)
(* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)
let m_Psi_tilde = Mat.init_cols n_ao n_ao (fun i j -> if j=q then m_Ksi_tilde.{i,2}
else if j=p then m_Ksi_tilde.{i,1}
else 0.);;
(* p q verif*)
(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *)
let m_Psi = Mat.init_cols n_ao n_ao (fun i j -> if j=q then m_Ksi.{i,2}
else if j=p then m_Ksi.{i,1}
else 0.);;
(* 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 m_interm = Mat.sub m_Phi m_Psi;;
(* Construction de la nouvelle matrice d'OMs phi~ *)
let m_Phi_tilde = Mat.add m_Psi_tilde m_interm;;
*)