25 KiB
25 KiB
None
<html>
<head>
</head>
</html>
In [ ]:
(* Algorithme de localisation des orbitales : Méthode de Boys ou Edminston-Ruedenberg *)
(*I : Position et base atomique*)
(* 1. Fonction création chaîne linéaire de n H *)
(* 2. Base atomique *)
(* II : Hartree-Fock *)
(* 1. Def pour HF *)
(* 2. Calcul de Hartree-Fock*)
(* III : Définition des fonctions pour la localisation *)
(* 1. Définitions de base nécessaire pour la suite *)
(* 2. Fonction de calcul de tous les alpha ER *)
(* 3. Fonction de calcul de tous les alpha Boys *)
(* 4. Pattern matching : calcul de alpha et de D *)
(* 5. Norme de alpha *)
(* 6. Détermination alpha_max et ses indices i et j *)
(* 7. Matrice de rotation 2 par 2 *)
(* 8. Fonction : i,j -> matrice Ksi (n par 2) *)
(* 9. Fonction : i~,j~ -> matrice Ksi~ (n par 2) *)
(* 10. Matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)
(* 11. Matrice intérmédiaire où les orbitales i et j ont été supprimées et remplacées par des 0*)
(* IV : Boucle de localisation des OMs *)
(* V : Fonctions pour la séparation et le rassemblement des matrices *)
(* 1. Fonction de création d'une list 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 *)
(* 4. Liste des OMs occupées *)
(* 5. Fonction de rassemblement de 2 matrices *)
(* 6. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)
(* VI : Calcul et assemblage des OMs occupées/virtuels *)
(* VII : Calcul Coef CI *)
(* 1. Def de la mo_basis pour le HF pre CI *)
(* 2. HF pour CI *)
(* 3. Calcul de coef CI *)
In [ ]:
open Lacaml.D
In [ ]:
(*I : Position et base atomique*)
(* 1. 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 20;;
In [ ]:
(* 2. Base atomique *)
let basis_string =
"
HYDROGEN
S 6
1 0.3552322122E+02 0.9163596281E-02
2 0.6513143725E+01 0.4936149294E-01
3 0.1822142904E+01 0.1685383049E+00
4 0.6259552659E+00 0.3705627997E+00
5 0.2430767471E+00 0.4164915298E+00
6 0.1001124280E+00 0.1303340841E+00
"
In [ ]:
(* II : Hartree-Fock *)
(* 1. Def pour HF *)
let nuclei =
Nuclei.of_xyz_string xyz_string
let basis =
Basis.of_nuclei_and_basis_string nuclei basis_string
let simulation =
Simulation.make ~charge:0 ~multiplicity:1 ~nuclei basis
let ao_basis =
Simulation.ao_basis simulation
let nocc =
let elec = Simulation.electrons simulation in
Electrons.n_alfa elec
In [ ]:
(* 2. Calcul de Hartree-Fock*)
let hf = HartreeFock.make simulation
In [ ]:
let mo_basis = MOBasis.of_hartree_fock hf
In [ ]:
let mo_coef = MOBasis.mo_coef mo_basis
In [ ]:
(* III : Définition des fonctions pour la localisation *)
(* 1. 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;;
In [ ]:
(* 2. Fonction de calcul de tous les alpha ER *)
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);;
In [ ]:
(* 3. Fonction de calcul de tous les 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.)));;
In [ ]:
(* 4. Pattern matching : calcul de alpha et de D *)
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;;
In [ ]:
(* 5. 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);;
In [ ]:
(* 6. 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 *)
type alphaij = {
alpha_max : float;
indice_ii : int;
indice_jj : int;};;
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
alpha_m.{i,j}
(*Printf.printf "%i %i\n%!" 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);;
In [ ]:
(* 7. 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 [ ]:
(* 8. 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} );;
In [ ]:
(* 9. 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;;
In [ ]:
(* 10. 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.)
In [ ]:
(* 11. 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
des coefs par la matrice *)
let f_interm m_C m_Psi = Mat.sub m_C m_Psi;;
In [ ]:
(* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)
let new_m_C m_C methode m_alpha epsilon =
(* 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
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 prev_m_alpha cc max_D=
(*Printf.printf "%i\n%!" n;
Printf.printf "%f\n%!" epsilon;*)
(*Util.debug_matrix "m_C" m_C;*)
(*Util.debug_matrix "new_alpha_m" prev_m_alpha;*)
if n == 0
then m_C
else
let m_new_m_C = new_m_C m_C methode prev_m_alpha epsilon
in
(* Fonction de pattern matching en fonction de la méthode *)
let alphad = m_alpha_d methode m_new_m_C
in
(* D critère à maximiser *)
let critere_D = alphad.d
in
let diff = max_D -. critere_D
in
let max_D =
if diff > 0.
then max_D
else critere_D
in
(*Printf.printf "%f\n%!" max_D;*)
Printf.printf "%e\n%!" critere_D;
(*Printf.printf "%f\n%!" diff;*)
(* Matrice des alphas *)
let m_alpha = alphad.m_alpha
in
(*let norme_alpha = norme m_alpha
in
Printf.printf "%f\n%!" norme_alpha;*)
if diff > 0. && epsilon >= 0.1
then localisation m_C methode (epsilon /. 2.) (n-1) critere_D prev_m_alpha cc max_D
else
let epsilon = 1.
in
(*Util.debug_matrix "new_alpha_m" m_alpha;*)
(*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 m_alpha cc max_D;;
In [ ]:
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
let max_D = 0.
in
localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc max_D;;
In [ ]:
(* V : Fonctions pour la séparation et le rassemblement des matrices *)
(* 1. 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;;
(* 2. 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;;
(* 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 *)
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);;
(* 4. 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;;
(* 5. Fonction de rassemblement de 2 matrices *)
let assemble m_occ m_vir = Mat.init_cols n_ao n_mo ( fun i j ->
if j > nocc
then m_vir.{i,j-(n_mo - nocc)}
else m_occ.{i,j});;
(* 6. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)
let m_occ , m_vir = split_mat m_C list_occ;;
In [ ]:
localise localisation m_occ "Boys" 1. 1000 10e-5;;