28 KiB
28 KiB
None
<html>
<head>
</head>
</html>
Sauvegarde des anciennes fonctions¶
In [ ]:
In [ ]:
In [ ]:
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 [ ]:
(*
(* 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;;
*)