10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-10-30 10:48:12 +01:00
QCaml/MOBasis/Localisation.ml

1199 lines
44 KiB
OCaml

(* Algorithme de localisation des orbitales *)
(* SOMMAIRE *)
(* 0. Paramètres du calcul *)
(* I. Définition des types *)
(* II : Définition des fonctions pour la localisation *)
(* 1. Définitions de base nécessaire pour la suite *)
(* 2. Fonction de calcul de la matrice des angles selon la méthode de Edminston *)
(* 3. Fonction calcul de la matrice des angles selon la méthode de Boys façon Edminston = NON FONCTIONNELLE *)
(* 4. Fonction de calcul de la matrice des angles selon la méthode de Boys *)
(* 5. Pattern matching : calcul de alpha et du critere D *)
(* 6. Norme d'une matrice *)
(* 7. Détermination de 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
Si alpha max < pi/2 on ajoute pi/2 à la matrice des alphas de manière récursive *)
(* 8. Matrice de rotation 2 par 2 *)
(* 9. Matrice de rotation n par n avec un angle fixe *)
(* 10. 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 *)
(* 11. 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*)
(* 12. Fonction pour la création de matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)
(* 13. 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 *)
(* III. Localisation des OMs *)
(* 1. Boucle de calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)
(* 2. Fonction de localisation *)
(* IV : Fonctions pour la séparation et le rassemblement des matrices *)
(* 1. Fonction de création d'une liste 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, 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 *)
(* 4. Liste des OMs occupées *)
(* 5. Liste des OMs de coeur *)
(* 6. Fonction de creation d'une liste des OMs pi virtuelles *)
(* 7. Fonction de rassemblement de 2 matrices *)
(* V. Calcul *)
(* 1. Séparation des OMs -> m_occ : OMs occupées; m_vir : OMs vacantes *)
(* 1.1 Liste des OMs pi virtuelles *)
(* 2. Application d'une légère rotation d'un angle fixé aux OMs pour améliorer la convergence *)
(* 3. Légende *)
(* 4. Rotation des orbitales : localise localisation new_X / "methode" / epsilon = 1. si <1 diminue l'angle des rotations / nb itération / critère de convergence *)
(* 4.1 Fonction de localisation de toutes les orbitales virtuelles sans separation *)
(* 4.2 Fonction de localisation en separant les orbitales pi et sigma *)
(* 4.3 Fonction de localisation en separant les orbitales pi et sigma et les 10 dernières sigma *)
(* 4.4 Fonction séparant les orbitales pi, la premiere et la seconde moitié des orbitales sigma *)
(* 5. Rassemblement des occupées et des virtuelles *)
(* 5.1 Rassemblement des orbitales occupées (coeur et non coeur) *)
(* 5.2 Pattern matching du type de separation des orbitales virtuelles et rassemblement des orbitales occupées et virtuelles *)
(* VI. Analyse *)
(* 1. Fonction de calcul de l'extension spatiale des orbitales *)
(* 2. Fonction de tri par ordre croissant d'une liste *)
(* 3. Fonction de tri par ordre croissant de l'extension spatiale des orbitales *)
(* 4. Moyenne/variance/ecart-type/mediane de l'extension spatiale des orbitales *)
(* 5. Fonction d'affichage *)
let methode = "Boys_er";; (* Method for th orbitals localization *)
let iteration = 1_000_000;; (* Maximum number of iteration for each localization *)
let cc = 10e-6;; (* Convergence criterion *)
let pre_angle = 0.001;; (* Rotation of the matrix with a small angle to improve the convergence *)
let epsilon = 1.;; (* Rotation angle = angle /. epsilon, default value : epsilon = 1. *)
let x_end = 1;; (* Choice of the x_end for the pi/sigma/end separation, size of the second matrix with the sigma orbitals *)
(* !!! WARNING : even if you don't use this localization method, 0 < x_end < total number of sigma orbitals.
If you don't use this method, use the default value, x_end = 1 WARNING !!! *)
(* Choice of the separation type of orbitals *)
(* For valence ones *)
type type_separation =
| Sigma_pi
| Occ_vir
| Custom of int list list
let string_of_separation = function
| Sigma_pi -> "Sigma/pi"
| Occ_vir -> "Occ/Vir"
| Custom _ -> "Custom"
open Lacaml.D
(* I. Types definitions *)
type alphaij = {
alpha_max : float;
indice_ii : int;
indice_jj : int;
}
type alphad = {
m_alpha : Mat.t;
d : float;
}
let localize mo_basis separation =
let simulation = MOBasis.simulation mo_basis in
let n_occ =
let elec = Simulation.electrons simulation
in
Electrons.n_alfa elec
in
let ao_basis = MOBasis.ao_basis mo_basis in
let n_core =
(*
(Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
*)
0
in
let charge = Simulation.charge simulation
|>Charge.to_int
in
(* II : Definition of the fonctions for the localization *)
(* 1. Basic definitions *)
let m_C = MOBasis.mo_coef mo_basis in
let n_ao = Mat.dim1 m_C in
let pi = acos(-1.) in
(* 2. Function for Edminston-Ruedenberg localization*)
(* Coefficient matrix n_ao x n_mo -> n_mo x n_mo matrix with angle(rad) between orbitals, Edminston localization criterion *)
(* matrix -> matrix, float *)
let f_alpha_er m_C =
let ee_ints = AOBasis.ee_ints ao_basis in
let n_mo = Mat.dim2 m_C 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);
(Mat.init_cols n_mo n_mo ( fun i j ->
if i= j
then 0.
else if asin( m_b12.{i,j} /. sqrt(m_b12.{i,j}**2. +. m_a12.{i,j}**2.)) > 0.
then 0.25 *. acos( -. m_a12.{i,j} /. sqrt(m_b12.{i,j}**2. +. m_a12.{i,j}**2.))
else -. 0.25 *. acos( -. m_a12.{i,j} /. sqrt(m_b12.{i,j}**2. +. m_a12.{i,j}**2.)))
,Vec.sum v_d)
in
(* 3. Function for the Boys localization like the Edminston-Ruedenberg localization *)
(* Coefficient matrix n_ao x n_mo -> n_mo x n_mo matrix with angle between orbitals, localization Boys like ER criterion *)
(* Matrix -> matrix, float *)
let f_alpha_boys_er m_C =
let multipoles = AOBasis.multipole ao_basis in
let n_mo = Mat.dim2 m_C in
let x =
Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let y =
Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let z =
Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let m_b12=
Mat.init_cols n_mo n_mo ( fun i j ->
(x.{i,i} -. x.{j,j}) *. x.{i,j} +.
(y.{i,i} -. y.{j,j}) *. y.{i,j} +.
(z.{i,i} -. z.{j,j}) *. z.{i,j})
in
let m_a12 =
Mat.init_cols n_mo n_mo ( fun i j ->
let x' = x.{i,i} -. x.{j,j}
and y' = y.{i,i} -. y.{j,j}
and z' = z.{i,i} -. z.{j,j} in
x.{i,j} *. x.{i,j} +.
y.{i,j} *. y.{i,j} +.
z.{i,j} *. z.{i,j} -.
0.25 *. (x'*.x' +. y'*.y' +. z'*.z')
)
in
(Mat.init_cols n_mo n_mo ( fun i j ->
if i=j
then 0.
else
let w = atan2 m_b12.{i,j} m_a12.{i,j} in
if w >= 0.
then 0.25 *. (pi -. w)
else -. 0.25 *. (pi +. w)
) , 0.5 *. (Vec.sum @@ Vec.init n_mo ( fun i ->
x.{i,i} *. x.{i,i} +.
y.{i,i} *. y.{i,i} +.
z.{i,i} *. z.{i,i} ) )
)
in
(* 4. Function for the original Boys localization *)
(* Coefficient matrix n_ao x n_mo -> (n_mo x n_mo matrix with angle between orbitals, localization criterion *)
(* Matrix -> matrix, float *)
let f_theta_boys m_C =
let multipoles = AOBasis.multipole ao_basis in
let n_mo = Mat.dim2 m_C
in
let x =
Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let y =
Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let z =
Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let r2 =
Multipole.matrix_r2 multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let m_b_0 =
Mat.init_cols n_mo n_mo (fun i j ->
r2.{i,i} +. r2.{j,j} -.
(x.{i,i} *. x.{i,i} +. x.{j,j} *. x.{j,j} +.
y.{i,i} *. y.{i,i} +. y.{j,j} *. y.{j,j} +.
z.{i,i} *. z.{i,i} +. z.{j,j} *. z.{j,j} )
)
in
let m_beta =
Mat.init_cols n_mo n_mo (fun i j ->
let x' = (x.{i,i} -. x.{j,j})
and y' = (y.{i,i} -. y.{j,j})
and z' = (z.{i,i} -. z.{j,j})
in
x'*.x' +. y'*.y' +. z'*.z' -.
4. *. ( x.{i,j} *. x.{i,j} +.
y.{i,j} *. y.{i,j} +.
z.{i,j} *. z.{i,j} )
)
in
let m_gamma =
Mat.init_cols n_mo n_mo (fun i j ->
4. *. ( x.{i,j} *. (x.{i,i} -. x.{j,j}) +.
y.{i,j} *. (y.{i,i} -. y.{j,j}) +.
z.{i,j} *. (z.{i,i} -. z.{j,j}) )
)
in
let m_theta = Mat.init_cols n_mo n_mo (fun i j ->
if i <> j then
0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j}
else 0.
)
in
let m_critere_B = Mat.init_cols n_mo n_mo (fun i j ->
0.5 *. (m_b_0.{i,j} +. 0.25 *. ((1. -. cos(4. *. m_theta.{i,j})) *. m_beta.{i,j} +. sin(4. *. m_theta.{i,j}) *. m_gamma.{i,j})))
in
m_theta, Vec.sum(Mat.as_vec m_critere_B)
in
(* Function to compute the new angle after rotation between the 2 orbitals j1 and j2, by the Boys like ER localization *)
(* n_mo x n_mo matrix with angle between orbitals, index of the orbital j1, index of the orbital j2 -> n_mo x 2 matrix with the new angles between j1 and j2, new localization criterion *)
(* matrix, integer, integer -> matrix, float *)
let f2_alpha_boys_er m_C j1 j2=
let multipoles = AOBasis.multipole ao_basis in
let n_mo = Mat.dim2 m_C
in
let x =
Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let y =
Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let z =
Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let m_b12 =
Mat.init_cols n_mo 2 (fun i j ->
if j = 1 then
(x.{i,i} -. x.{j1,j1}) *. x.{i,j1} +.
(y.{i,i} -. y.{j1,j1}) *. y.{i,j1} +.
(z.{i,i} -. z.{j1,j1}) *. z.{i,j1}
else
(x.{i,i} -. x.{j2,j2}) *. x.{i,j2} +.
(y.{i,i} -. y.{j2,j2}) *. y.{i,j2} +.
(z.{i,i} -. z.{j2,j2}) *. z.{i,j2}
)
in
let m_a12 =
Mat.init_cols n_mo 2 (fun i j ->
if j = 1 then
let x' = (x.{i,i} -. x.{j1,j1})
and y' = (y.{i,i} -. y.{j1,j1})
and z' = (z.{i,i} -. z.{j1,j1}) in
x.{i,j1} *. x.{i,j1} +.
y.{i,j1} *. y.{i,j1} +.
z.{i,j1} *. z.{i,j1} -.
0.25 *. (x' *. x' +. y' *. y' +. z' *. z' )
else
let x' = (x.{i,i} -. x.{j2,j2})
and y' = (y.{i,i} -. y.{j2,j2})
and z' = (z.{i,i} -. z.{j2,j2}) in
x.{i,j2} *. x.{i,j2} +.
y.{i,j2} *. y.{i,j2} +.
z.{i,j2} *. z.{i,j2} -.
0.25 *. (x' *. x' +. y' *. y' +. z' *. z' )
)
in
(Mat.init_cols n_mo 2 ( fun i j ->
if i=j1 && j=1
then 0.
else if i = j2 && j=2
then 0.
else
let x = atan2 m_b12.{i,j} m_a12.{i,j} in
if x >= 0. then 0.25 *. (pi -. x)
else -. 0.25 *. (pi +. x)
),
Vec.sum @@ Vec.init n_mo ( fun i ->
x.{i,i} *. x.{i,i} +.
y.{i,i} *. y.{i,i} +.
z.{i,i} *. z.{i,i} )
)
in
(* Function to compute the new angle after rotation between the 2 orbitals j1 and j2, by the Boys localization *)
(* n_mo x n_mo matrix with angles between orbitals, index of the orbital j1, index of the orbital j2 -> n_mo x 2 matrix with the new angles between j1 and j2, new localization criterion *)
(* matrix, integer, integer -> matrix, float *)
let f2_theta_boys m_C j1 j2 =
let multipoles = AOBasis.multipole ao_basis in
let n_mo = Mat.dim2 m_C
in
let x =
Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let y =
Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let z =
Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let r2 =
Multipole.matrix_r2 multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in
let m_beta =
Mat.init_cols n_mo 2 (fun i j ->
if j = 1 then
let x' = (x.{i,i} -. x.{j1,j1})
and y' = (y.{i,i} -. y.{j1,j1})
and z' = (z.{i,i} -. z.{j1,j1}) in
x' *. x' +. y' *. y' +. z' *. z'
-. 4. *. ( x.{i,j1} *. x.{i,j1} +.
y.{i,j1} *. y.{i,j1} +.
z.{i,j1} *. z.{i,j1} )
else
let x' = (x.{i,i} -. x.{j2,j2})
and y' = (y.{i,i} -. y.{j2,j2})
and z' = (z.{i,i} -. z.{j2,j2}) in
x' *. x' +. y' *. y' +. z' *. z'
-. 4. *. ( x.{i,j2} *. x.{i,j2} +.
y.{i,j2} *. y.{i,j2} +.
z.{i,j2} *. z.{i,j2} )
)
in
let m_gamma =
let gamma g = Mat.init_cols n_mo 2 (fun i j ->
if j = 1
then 4. *. g.{i,j1} *. (g.{i,i} -. g.{j1,j1})
else 4. *. g.{i,j2} *. (g.{i,i} -. g.{j2,j2}))
in Mat.add (gamma x) (Mat.add (gamma y) (gamma z))
in
let m_theta = Mat.init_cols n_mo 2 (fun i j ->
if i=j1 && j=1
then 0.
else if i = j2 && j=2
then 0.
else +. 0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j})
in
let m_critere_B = Vec.init n_mo (fun i ->
r2.{i,i} -. x.{i,i}*.x.{i,i} -. y.{i,i}*.y.{i,i} -. z.{i,i}*.z.{i,i})
in m_theta, Vec.sum( m_critere_B)
in
(* 5. Pattern matching : choice of the localization method *)
(* The method, coefficient matrix n_ao x n_mo -> n_mo x n_mo matrix with the angle between orbitals using previous function, localization criterion *)
(* string, matrix -> matrix, float *)
let m_alpha_d methode m_C =
let alpha methode =
match methode with
| "Boys_er"
| "BOYS_ER"
| "boys_er" -> let alpha_boys , d_boys = f_alpha_boys_er m_C in {m_alpha = alpha_boys; d = d_boys}
| "Boys"
| "boys"
| "BOYS" -> let theta_boys, b_boys = f_theta_boys m_C in {m_alpha = theta_boys; d = b_boys}
| "ER"
| "er" -> let alpha_er , d_er = f_alpha_er m_C in {m_alpha = alpha_er; d = d_er}
| _ -> invalid_arg "Unknown method, please enter Boys or ER"
in
alpha methode
in
(* 5. Pattern matching : choice of the localization method to compute the n_mo x 2 new angles and not the n_mo x n_mo angles, after the rotation of i and j orbitals *)
(* The method, coefficient matrix n_ao x n_mo, index of the orbital i , index of the orbital j -> n_mo x 2 matrix with the new angle between orbitals i and j using previous function, new localization criterion *)
(* string, matrix, integer, integer -> matrix, float *)
let f_speed_alpha methode m_C indice_i indice_j =
let alpha_speed methode =
match methode with
| "Boys_er"
| "boys_er" -> let m_alpha, critere_D = f2_alpha_boys_er m_C indice_i indice_j in {m_alpha = m_alpha; d = critere_D}
| "BOYS" -> let m_alpha, critere_D = f2_theta_boys m_C indice_i indice_j in {m_alpha = m_alpha; d = critere_D}
| _ -> invalid_arg "Unknown method, please enter Boys or ER"
in
alpha_speed methode
in
(* 6. Norm of a matrix *)
(* Matrix -> float *)
let norm m =
sqrt (Mat.syrk_trace m)
in
(* 7. Calculation of the max angle[-pi/.4 , pi/.4] and its indexes i and j in the matrix *)
(* n_mo x n_mo matrix with angles between orbitals, coefficient matrix n_ao x n_mo, max number of iteration -> max angle, line location, column location *)
(* matrix, matrix, integer -> float, integer, integer *)
let rec new_m_alpha m_alpha m_C n_rec_alpha=
let pi_2 = 0.5 *. pi in
let n_mo = Mat.dim2 m_C in
let alpha_m =
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 m_alpha.{i,j}
)
in
(* Location of the max angle in the matrix *)
(* angle matrix -> location of the max element *)
(* matrix -> integer *)
let max_element alpha_m =
Mat.as_vec alpha_m
|> iamax
in
(* Indexes i and j of the max angle *)
(* -> integer, integer *)
let indice_ii, indice_jj =
let max = max_element alpha_m
in
(max - 1) mod n_mo +1, (max - 1) / n_mo +1
in
(* Value of the max angle *)
(* angle matrix -> value of the max angle*)
(* matrix -> float *)
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
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
(* 8. Rotation matrix 2x2 *)
(* angle -> 2x2 rotation matrix *)
(* float -> matrix *)
let f_R alpha =
[| [| cos alpha ; -. sin alpha |] ;
[| sin alpha ; cos alpha |] |]
|> Mat.of_array
in
(* 9. Rotation matrix nxn *)
(* Coefficient matrix n_ao x n_mo, angle -> new coefficient matrix n_ao x n_mo after a rotation with a small angle *)
(* matrix, float -> matrix *)
let rotate m_C angle=
let n_mo = Mat.dim2 m_C in
(* Antisymmetrization of the matrix and diagonal terms = 0, X *)
(* angle -> Antisymmetric matrix with this angle and diagonal terms = 0 *)
(* float -> matrix *)
let m_angle angle=
Mat.init_cols n_mo n_mo (fun i j ->
if i = j
then 0.
else if i>j
then -. angle
else angle)
in
let m = m_angle angle in
(* Square of the matrix X² *)
(* -> X . X *)
(* -> matrix*)
let mm = gemm m m in
(* Diagonalization of X² *)
(* diag = vector that contains the eigenvalues (-tau²), mm contains the eigenvectors (W) *)
(* -> vector *)
let diag = syev ~vectors:true mm in
(* -tau² to tau² *)
(* -> vector that contains the eigenvalues (tau²) *)
(* -> vector *)
let square_tau= Vec.abs diag in
(* tau² to tau *)
(* -> vector that contains the eigenvalues tau *)
(* -> vector *)
let tau = Vec.sqrt square_tau in
(* Calculation of cos tau from the vector tau *)
(* -> cos(tau) *)
(* -> float *)
let cos_tau = Vec.cos tau in
(* Matrix cos tau *)
(* -> matrix with diagonal terms cos(tau) and 0 also *)
(* -> matrix *)
let m_cos_tau =
Mat.init_cols n_mo n_mo (fun i j ->
if i=j
then cos_tau.{i}
else 0.)
in
(* Calcul of sin(tau) from the vector tau *)
let sin_tau = Vec.sin tau in
(* Matrix sin tau *)
(* -> matrix with diagonal terms sin(tau) and 0 also *)
(* -> matrix *)
let m_sin_tau =
Mat.init_cols n_mo n_mo (fun i j ->
if i=j
then sin_tau.{i}
else 0.)
in
(* Calculation of the vector tau^-1 *)
(* -> vector tau^-1 *)
(* -> vector *)
let tau_1 =
Vec.init n_mo (fun i ->
if tau.{i}<= 0.001
then 1.
else 1. /. tau.{i})
in
(* Calculation of the matrix tau^⁻1 *)
(* -> matrix tau^-1 *)
(* -> matrix *)
let m_tau_1 =
Mat.init_cols n_mo n_mo (fun i j ->
if i=j
then tau_1.{i}
else 0.)
in
(* W cos(tau) Wt *)
(* -> matrix *)
let a = gemm mm @@ gemm ~transb:`T m_cos_tau mm in
(* W tau^-1 sin(tau) Wt X *)
(* -> matrix *)
let b = gemm mm @@ gemm m_tau_1 @@ gemm m_sin_tau @@ gemm ~transa:`T mm m in
(* Sum of a + b -> Rotation matrix *)
(* -> Rotation matrix *)
(* -> matrix *)
let m_r = Mat.add a b
in
gemm m_C m_r
in
(* Function to extract 2 vectors i and j in a matrix *)
(* Coefficient matrix n_ao x n_mo, index of the orbital i, index of the orbital j -> n_mo x 2 matrix with the eigenvetors i and j*)
(* matrix, integer, integer -> matrix *)
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
(* Function to apply a rotation with a 2x2 rotation matrix on a n_mo x 2 matrix that contains i and j *)
(* 2x2 rotation matrix, n_mo x 2 eigenvetors(i and j) matrix -> n_mo x 2 new eigenvetors(i~ and j~) after rotation *)
(* matrix, matrix -> matrix *)
let f_Ksi_tilde m_R m_Ksi = gemm m_Ksi m_R
in
(* Function to create intermediate coefficient matrix in order to delete the i j orbitals and add i~ j~ in the coefficient matrix *)
(* n_mo x 2 matrix with k and l orbitals, index of k, index of l, n_ao x n_mo matrix coefficient -> matrix where all terms are equal to 0 except the k and l columns that contain the eigenvetors k and l *)
(* matrix, integer, integer, matrix -> matrix *)
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
(* Function to create a intermediate angle matrix in order to delete all the previous angle between i and all the other orbitals, j and all the other orbitals.
And add the new angles between these orbitals *)
(* n_mo x 2 matrix with angles between k and all the other orbitals l and all the other orbitals, index of k, index of l, coefficent matrix n_ao x n_mo -> matrix where all terms are equal to 0 except the k and l lines and columns that contains the angles between the orbitals *)
(* matrix, integer, integer, matrix -> matrix *)
let f_k_angle mat indice_i indice_j m_C =
let n_mo = Mat.dim2 m_C
in
Mat.init_cols n_mo n_mo (fun i j ->
if j=indice_i
then mat.{i,1}
else if j=indice_j
then mat.{i,2}
else if i = indice_i
then -. mat.{j,1}
else if i = indice_j
then -. mat.{j,2}
else 0.)
in
(* Intermediate matrix where the i and j orbitals are equal to 0 *)
(* Coefficient matrix n_ao x n_mo, matrix where terms are equal to 0 except the k and l columns that contain the eigenvetors k and l *)
(* matrix, matrix -> matrix *)
let f_interm m_C mat =
Mat.sub m_C mat
in
(* Function to compute the new coefficient matrix after a rotation between 2 orbitals *)
(* coefficient matrix n_ao x n_mo, angle matrix, parameter -> new coefficient matrix n_ao x n_mo *)
let new_m_C m_C m_alpha epsilon =
(* alphaij contains the max angle of the angle matrix and the indexes *)
let n_rec_alpha = 10 in
let alphaij = new_m_alpha m_alpha m_C n_rec_alpha in
(* Value of the angle *)
let alpha = (alphaij.alpha_max) *. epsilon in
(* Localtion of the max angle *)
let indice_i = alphaij.indice_ii in
let indice_j = alphaij.indice_jj in
(* Rotation matrix *)
let m_R = f_R alpha in
(* Matrix with i and j orbitals *)
let m_Ksi = f_Ksi indice_i indice_j m_C in
(* Matrix with the new i~ and j~ orbitals *)
let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi in
(* Matrix to delete the i and j orbitals in the coefficient matrix *)
let m_Psi = f_k m_Ksi indice_i indice_j m_C in
(* Matrix to add the i~ and j~ orbitals in the coefficient matrix *)
let m_Psi_tilde = f_k m_Ksi_tilde indice_i indice_j m_C in
(* Coefficient matrix without the i and j orbitals *)
let m_interm = f_interm m_C m_Psi
in
(* Matrice after rotation, max angle, index i, index j *)
(Mat.add m_Psi_tilde m_interm, alpha, indice_i, indice_j)
in
(* Function to compute the new angle matrix after one rotation between 2 orbitals *)
(* Previous angle matrix n_mo x n_mo, new coefficient matrix n_ao x n_mo, previous angle matrix, index of the orbital i, index of the orbital j -> new angle matrix, new localization criterion *)
(* matrix, matrix, matrix, integer, integer -> matrix, float *)
let m_alpha m_new_m_C prev_m_alpha indice_i indice_j =
(* n_mo x 2 matrix that contains the new angles between i,j with the other orbitals *)
let speed_alphad = f_speed_alpha methode m_new_m_C indice_i indice_j in
(* New localization criterion *)
let critere_D = speed_alphad.d in
(* Previous angles *)
let ksi_angle = f_Ksi indice_i indice_j prev_m_alpha in
(* New angles *)
let ksi_angle_tilde = speed_alphad.m_alpha in
(* Matrix to delete the previous angles *)
let psi_angle = f_k_angle ksi_angle indice_i indice_j m_new_m_C in
(* Matrix to add the new angles *)
let psi_tilde_angle = f_k_angle ksi_angle_tilde indice_i indice_j m_new_m_C in
(* Matrix without the angles between i, j and the orbitals *)
let m_interm = f_interm prev_m_alpha psi_angle
in
(* New angle matrix, new localization criterion *)
(Mat.add m_interm psi_tilde_angle, critere_D)
in
(* III. Localization *)
(* Loop to compute the new coefficient matrix avec n orbital rotations *)
(* Coefficeient matrix n_ao x n_mo, localization method, parameter, max number of iteration, previous localization criterion, previous angle matrix, convergence criterion -> new coefficient matrix after n orbital rotations *)
(* matrix, string, float, integer, float, matrix, float -> matrix *)
let rec localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc =
if n == 0
then m_C
else
(* New coefficient matrix after one rotation *)
let m_new_m_C, alpha_max, indice_i, indice_j = new_m_C m_C prev_m_alpha epsilon in
(* New angle matrix after one rotation *)
let m_alpha, critere_D = m_alpha m_new_m_C prev_m_alpha indice_i indice_j in
(* Matrix norm, angle average *)
let norm_alpha = norm m_alpha in
let moyenne_alpha = norm_alpha /. (float_of_int((Mat.dim1 m_C)* (Mat.dim2 m_C))) in
let alpha = alpha_max /. epsilon in
Printf.printf "%i %f %f %f %f\n%!" (iteration-n) critere_D alpha norm_alpha moyenne_alpha;
(* Convergence criterion *)
if alpha_max**2. < cc**2.
then m_new_m_C
else localisation m_new_m_C methode epsilon (n-1) critere_D m_alpha cc
in
(* 2. Localization function *)
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
localisation m_C methode epsilon n prev_critere_D prev_m_alpha cc
in
(* IV : Functions to split and assemble matrix *)
(* Function to create a list from the missing elements of an other list*)
(* Coefficient matrix n_ao x n_mo, list -> list of missing element of the previous list *)
(* matrix, integer list -> integer list *)
let miss_elem mat lst =
let n_mo = Mat.dim2 mat in
Util.list_range 1 n_mo
|> List.map (fun i ->
if List.mem i lst
then 0
else i)
|> List.filter (fun x -> x > 0)
in
(* Function to split a matrix in 2 matrix, the first matrix corresponds to the column number whose number is in the list, the second matrix corresponds to the column which are not in the list *)
(* Coefficient matrix n_ao x n_mo, list -> matrix, matrix *)
(* matrix, integer list -> matrix, matrix*)
let split_mat mat lst =
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 lst in
let list_2 = miss_elem mat lst 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)
in
(* Function to create a list of pi MOs *)
(* Coefficient matrix n_ao x n_mo -> list of pi orbitals *)
(* matrix -> integer list *)
let list_pi mat =
let n_mo = Mat.dim2 mat in
let m_pi = Mat.init_cols n_ao (n_mo+1) ( fun i j ->
if j=1 && abs_float m_C.{i,j} < 1.e-8
then 0.
else if j=1
then abs_float m_C.{i,j}
else if abs_float mat.{i,j-1} < 1.e-8
then 0.
else abs_float mat.{i,j-1})
in
let vec = Mat.to_col_vecs m_pi in
(* TODO gemv *)
let vec_dot = Vec.init (Mat.dim2 mat) (fun i -> dot vec.(0) vec.(i)) in
let list_pi =
Util.list_range 1 (Mat.dim2 mat)
|> List.map (fun i -> if abs_float vec_dot.{i} = 0. then i else 0)
in
let rec remove x = function
| [] -> []
| var :: tail -> if var = x
then remove x tail
else var :: (remove x tail)
in
remove 0 list_pi
in
(* Function to assemble to matrix with the same number of line *)
(* Coefficient matrix n_ao x n, coefficient matrix n_ao x m -> Coefficient matrix n_ao x (n+m) *)
(* matrix, matrix -> matrix*)
let assemble mat_1 mat_2 =
let v_mat_1 = Mat.to_col_vecs mat_1 in
let v_mat_2 = Mat.to_col_vecs mat_2 in
let m = Array.append v_mat_1 v_mat_2
in
Mat.of_col_vecs m
in
(* V. Calculation *)
(* 1. Caption *)
let text = "Method : " ^ methode ^ "\nSeparation type : "
^ (string_of_separation separation)
^ "\nMax number of iterations : " ^ string_of_int(iteration) ^ "\ncc : "
^ string_of_float(cc)^"\nepsilon : "^string_of_float(epsilon)^"\ncharge : "
^ string_of_int(charge)^"\nRotation pre angle : "^string_of_float(pre_angle)
in
let caption = "Criterion max_angle norm_angle average_angle"
in
Printf.printf "%s\n" "";
Printf.printf "%s\n" "Orbital localization";
Printf.printf "%s\n" "";
Printf.printf "%s\n" text;
Printf.printf "%s\n" "";
Printf.printf "%s\n" caption;
Printf.printf "%s\n" "";
(* Localization function *)
(* Coefficient matrix n_ao x n_mo -> localized coefficient matrix n_ao x n_mo *)
(* matrix -> matrix *)
let f_localise mat =
if Mat.dim2 mat = 0
then Printf.printf "%s\n" "0 orbital, no localization";
let new_mat =
if Mat.dim2 mat = 0
then mat
else rotate mat pre_angle
in
if Mat.dim2 mat = 0
then new_mat
else localise localisation new_mat methode epsilon iteration cc
in
(* Localization function with n separations expressed as a list of lists of
orbital indices *)
let f_loc_list mat list_of_lists =
let vec_list = Mat.to_col_vecs mat in
list_of_lists
|> List.map (fun lst ->
lst
|> List.map (fun col -> vec_list.(col-1))
|> Mat.of_col_vecs_list
|> f_localise
|> Mat.to_col_vecs_list
)
|> List.concat
|> Mat.of_col_vecs_list
in
let n_mo = Mat.dim2 m_C in
let list_full = Util.list_range 1 n_mo in
let list_core = Util.list_range 1 n_core in
let list_pi = list_pi m_C in
let list_sigma = List.filter (fun i -> not (List.mem i list_pi)) list_full in
let list_pi_occ = List.filter (fun i -> i > n_core && i <= n_occ) list_pi in
let list_sigma_occ = List.filter (fun i -> i > n_core && i <= n_occ) list_sigma in
let list_pi_vir = List.filter (fun i -> i > n_occ) list_pi in
let list_sigma_vir = List.filter (fun i -> i > n_occ) list_sigma in
let list_occ = Util.list_range (n_core+1) n_occ in
let list_vir = Util.list_range (n_occ+1) n_mo in
Printf.printf "Orbital lists:";
Printf.printf "\nCore: [";
List.iter (fun i -> Printf.printf "%d," i) list_core;
Printf.printf "\b]\nOccupied: [";
List.iter (fun i -> Printf.printf "%d," i) list_occ;
Printf.printf "\b]\nVirtual: [";
List.iter (fun i -> Printf.printf "%d," i) list_vir;
Printf.printf "\b]\nSigma occ: [";
List.iter (fun i -> Printf.printf "%d," i) list_sigma_occ;
Printf.printf "\b]\nPi occ: [";
List.iter (fun i -> Printf.printf "%d," i) list_pi_occ;
Printf.printf "\b]\nSigma vir: [";
List.iter (fun i -> Printf.printf "%d," i) list_sigma_vir;
Printf.printf "\b]\nPi vir: [";
List.iter (fun i -> Printf.printf "%d," i) list_pi_vir;
Printf.printf "\b]\n\n";
(* Separation of the occupied and virtual MOs, valence and core MOs *)
let m_occ , m_vir = split_mat m_C list_occ in
let m_core, m_occu = split_mat m_occ list_core in
let m_assemble_loc =
let list_of_lists =
match separation with
| Occ_vir -> [ list_core ; list_occ ; list_vir ]
| Sigma_pi -> [ list_core ; list_sigma_occ ; list_pi_occ ; list_pi_vir ; list_sigma_vir ]
| Custom l -> l
in
f_loc_list m_C list_of_lists
in
let m_assemble_occ, m_assemble_vir = split_mat m_assemble_loc list_occ in
let m_assemble_core, m_assemble_occu = split_mat m_assemble_occ list_core in
(* VI. Analysis *)
(* Function to compute the spatial extent of MOs *)
(* Coefficient matrix n_ao x n_mo -> vector where the i th element is equal to the spatial extent of the i th orbital *)
(* Matrix -> vector *)
let distrib mat=
if Mat.dim1 mat = 0
then Mat.as_vec mat
else
(let multipoles = AOBasis.multipole ao_basis in
let n_mo = Mat.dim2 mat in
let x =
Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
in
let y =
Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
in
let z =
Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
in
let r2 =
Multipole.matrix_r2 multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat
in
Vec.init n_mo (fun i ->
r2.{i,i} -. x.{i,i}*.x.{i,i} -. y.{i,i}*.y.{i,i} -. z.{i,i}*.z.{i,i}))
in
(* Sorting function *)
(* tri insertion (<) list -> list in ascending order *)
(* parameter, list -> list *)
let rec insere comparaison elem liste = match liste with
| [] -> elem::[]
| tete::queue -> if comparaison elem tete
then elem :: liste
else tete :: insere comparaison elem queue
in
let rec tri_insertion comp = function
| [] -> []
| tete::queue -> insere comp tete (tri_insertion comp queue)
in
(* Sorting function to sort orbital spatial extent by ascending order *)
(* coefficient matrix n_ao x n_mo -> matrix n_ao x 2, first column corresponds to the spatial extent in function of the MOs number, the second column corresponds to the spatial extent in ascending order *)
(* matrix -> matrix *)
let m_distribution mat =
let vec_distrib= distrib mat in
let list_tri_distrib = tri_insertion (<) (Vec.to_list vec_distrib) in
let vec_tri_distrib = Vec.of_list list_tri_distrib
in
Mat.init_cols (Vec.dim vec_distrib) 2 (fun i j ->
if j=1
then vec_distrib.{i}
else vec_tri_distrib.{i})
in
(* Average/variance/standart deviation/median of spatial extent *)
(* Function to compute the average *)
let f_average mat =
if Mat.dim1 mat = 0
then 0.
else Vec.sum (distrib mat) /. float_of_int(Vec.dim (distrib mat))
in
(* Function to compute the variance *)
let f_variance mat =
if (Vec.dim (distrib mat)) <= 1
then 0.
else Vec.sum(Vec.init (Vec.dim (distrib mat)) (fun i ->
((distrib mat).{i}-. (f_average mat))**2.)) /. float_of_int(Vec.dim (distrib mat) - 1)
in
(* Function to compute the standard deviation *)
let f_stand_dev mat =
if (Vec.dim (distrib mat)) <= 1
then 0.
else sqrt(abs_float(f_variance mat))
in
(* Fonction de calcul de la mediane *)
let f_median mat =
if (Vec.dim (distrib mat)) = 0
then 0.
else
let vec_distrib= distrib mat in
let list_tri_distrib = tri_insertion (<) (Vec.to_list vec_distrib) in
let vec_tri_distrib = Vec.of_list list_tri_distrib
in
if (Vec.dim vec_distrib) mod 2 = 0
then (vec_tri_distrib.{(Vec.dim vec_distrib)/2} +. vec_tri_distrib.{((Vec.dim vec_distrib)/2)+1}) /. 2.
else vec_tri_distrib.{int_of_float(float_of_int(Vec.dim vec_distrib)/. 2.) + 1}
in
(* Display function *)
(* Coefficient matrix n_ao x n_mo -> string *)
(* matrix -> string *)
let analyse mat =
let average = f_average mat in
let variance = f_variance mat in
let stand_dev = f_stand_dev mat in
let median = f_median mat
in
"Average : "^string_of_float(average)^"; Median : "^string_of_float(median)^"; Variance : "^string_of_float(variance)^"; Standard deviation : "^string_of_float(stand_dev)
in
(* Display the numer of MOs *)
(* Coefficient matrix n_ao x n_mo -> string *)
(* matrix -> string *)
let n_orb mat = "Number of molecular orbitals : "^string_of_int(Mat.dim2 mat) in
(* Display *)
Printf.printf "%s\n" "";
Printf.printf "%s\n" (n_orb m_C);
Printf.printf "%s %i\n" "Number of pi orbitals : " (List.length list_pi );
Printf.printf "%s %i\n" "Number of occupied pi orbitals : " (List.length list_pi_occ );
Printf.printf "%s %i\n" "Number of virtual pi orbitals : " (List.length list_pi_vir );
Printf.printf "%s %i\n" "Number of sigma orbitals : " (List.length list_sigma);
Printf.printf "%s %i\n" "Number of occupied sigma orbitals : " (List.length list_sigma_occ);
Printf.printf "%s %i\n" "Number of virtual sigma orbitals : " (List.length list_sigma_vir);
Printf.printf "%s\n" "";
Printf.printf "%s %s\n" "All MOs before the localization" (analyse m_C);
Printf.printf "%s %s\n" "All MOs after the localization" (analyse m_assemble_loc);
Printf.printf "%s\n" "";
Util.debug_matrix "Distribution of the spatial extent, before localization / sorting before localization / after localization / sorting after localization"
(assemble(m_distribution m_C) (m_distribution m_assemble_loc));
Printf.printf "%s\n" "";
Printf.printf "%s %s\n" "Occupied orbitals : " (n_orb m_occ);
Printf.printf "%s %s\n" "Occupied orbitals before localization" (analyse m_occ);
Printf.printf "%s %s\n" "Occupied orbitals after localization" (analyse m_assemble_occ);
Printf.printf "%s\n" "";
Printf.printf "%s\n" "";
Printf.printf "%s %s\n" "Core orbitals : " (n_orb m_core);
Printf.printf "%s %s\n" "Core orbitalses before localization : " (analyse m_core);
Printf.printf "%s %s\n" "Core orbitals after localization : " (analyse m_assemble_core);
Printf.printf "%s\n" "";
Printf.printf "%s\n" "";
Printf.printf "%s %s\n" "Valence orbitals : " (n_orb m_occu);
Printf.printf "%s %s\n" "Valence orbitals before localization : " (analyse m_occu);
Printf.printf "%s %s\n" "Valence rbitals after localization : " (analyse m_assemble_occu);
Printf.printf "%s\n" "";
Printf.printf "%s\n" "";
Printf.printf "%s %s\n" "Virtual orbitals : " (n_orb m_vir);
Printf.printf "%s %s\n" "Virtual orbitals before localization : " (analyse m_vir);
Printf.printf "%s %s\n" "Virtual orbitals after localization : " (analyse m_assemble_vir);
Printf.printf "%s\n" "";
m_assemble_loc