From 158d149036e1ec3beaa34ad884903b1b3c352013 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 9 Jul 2020 23:25:38 +0200 Subject: [PATCH] Optimizations in localization --- Basis/Multipole.ml | 5 +- Basis/Multipole.mli | 3 + MOBasis/Localisation.ml | 454 ++++++++++++++++++++-------------------- 3 files changed, 231 insertions(+), 231 deletions(-) diff --git a/Basis/Multipole.ml b/Basis/Multipole.ml index 2491902..9cf5aa0 100644 --- a/Basis/Multipole.ml +++ b/Basis/Multipole.ml @@ -21,6 +21,7 @@ let matrix_z t = t.(2) let matrix_x2 t = t.(3) let matrix_y2 t = t.(4) let matrix_z2 t = t.(5) +let matrix_r2 t = Mat.add t.(3) @@ Mat.add t.(4) t.(5) let matrix_xy t = t.(6) let matrix_yz t = t.(7) let matrix_zx t = t.(8) @@ -215,8 +216,8 @@ let of_basis basis = done; done; done; - for k=0 to 14 do - Mat.detri result.(k); + for k=0 to Array.length result - 1 do + Mat.detri result.(k) done; result diff --git a/Basis/Multipole.mli b/Basis/Multipole.mli index 74346bd..43564fa 100644 --- a/Basis/Multipole.mli +++ b/Basis/Multipole.mli @@ -31,6 +31,9 @@ val matrix_y2 : t -> Mat.t val matrix_z2 : t -> Mat.t (** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *) +val matrix_r2 : t -> Mat.t +(** {% $$ \langle \chi_i | r^2 | \chi_j \rangle $$ %} *) + val matrix_x3 : t -> Mat.t (** {% $$ \langle \chi_i | x^3 | \chi_j \rangle $$ %} *) diff --git a/MOBasis/Localisation.ml b/MOBasis/Localisation.ml index e7d1148..de7a95a 100644 --- a/MOBasis/Localisation.ml +++ b/MOBasis/Localisation.ml @@ -80,7 +80,7 @@ let x_end = 1;; (* Choice of the x_end for the pi/sigma/end separation, size of (* For valence ones *) let type_separation_occu = "pi/sigma";; (* Separation of the pi and sigma orbitals *) -(*let type_separation_occu = "1/2";;*) (* Separation of the first and second half of the orbitals *) +(*let type_separation_occu = "1/2";;*) (* Separation of the first and second half of the orbitals *) (*let type_separation_occu = "pi/sigma/end";;*) (* Separation pi,(n-x_end) sigma and x_end sigma *) (*let type_separation_occu = "pi/1/2";;*) (* Separation pi, (n/2) first sigma and (n/2) last sigma *) (*let type_separation_occu = "full";;*) (* Without separation *) @@ -120,11 +120,11 @@ let localize mo_basis = (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 in let charge = Simulation.charge simulation - |>Charge.to_int + |>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 @@ -134,14 +134,14 @@ let localize mo_basis = (* 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) @@ -149,7 +149,7 @@ let localize mo_basis = 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 -> @@ -159,44 +159,44 @@ let localize mo_basis = ) (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) = = \sum_p

C_{pi} *) gemm ~transa:`T ~c:m_qr_i m_p_qr m_C in - + let m_q_ri = (* Transformation de la matrice (qr,i) en (q,ri) *) Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo) in - + let m_ri_j = (* (ri,j) = = \sum_q C_{bj} *) gemm ~transa:`T ~c:m_ri_j m_q_ri m_C in - + let m_r_ij = (* Transformation de la matrice (ri,j) en (r,ij) *) Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo) in - + let m_ij_k = (* (ij,k) = = \sum_r C_{rk} *) gemm ~transa:`T ~c:m_ij_k m_r_ij m_C in - + let m_ijk = (* Transformation de la matrice (ei,j) en (e,ij) *) Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |] |> Bigarray.array3_of_genarray in - + Array.iter (fun j -> Array.iter (fun i -> m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j}); @@ -207,51 +207,51 @@ let localize mo_basis = 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 + 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 f_alpha_boys_er m_C = + let multipoles = AOBasis.multipole ao_basis in let n_mo = Mat.dim2 m_C in let phi_x_phi = - Multipole.matrix_x multipoles + Multipole.matrix_x multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C - in + in let phi_y_phi = - Multipole.matrix_y multipoles + Multipole.matrix_y multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in let phi_z_phi = - Multipole.matrix_z multipoles + Multipole.matrix_z multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C - in - let m_b12= + 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 + in Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi)) - in + in let m_a12 = - let a12 g = Mat.init_cols n_mo n_mo ( fun i j -> + 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 + (Mat.init_cols n_mo n_mo ( fun i j -> + if i=j then 0. else if +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} >= 0. then pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} @@ -260,65 +260,65 @@ let localize mo_basis = 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 = +(* 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 phi_x_phi = - Multipole.matrix_x multipoles + Multipole.matrix_x multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C - in + in let phi_y_phi = - Multipole.matrix_y multipoles + Multipole.matrix_y multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in let phi_z_phi = - Multipole.matrix_z multipoles + Multipole.matrix_z multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in let phi_x2_phi = - Multipole.matrix_x2 multipoles + Multipole.matrix_x2 multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C - in + in let phi_y2_phi = - Multipole.matrix_y2 multipoles + Multipole.matrix_y2 multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in let phi_z2_phi = - Multipole.matrix_z2 multipoles + Multipole.matrix_z2 multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in let m_b_0 = let b_0 g h = Mat.init_cols n_mo n_mo (fun i j -> h.{i,i} +. h.{j,j} -. (g.{i,i}**2. +. g.{j,j}**2.)) in - Mat.add (b_0 phi_x_phi phi_x2_phi) (Mat.add (b_0 phi_y_phi phi_y2_phi) (b_0 phi_z_phi phi_z2_phi)) + Mat.add (b_0 phi_x_phi phi_x2_phi) (Mat.add (b_0 phi_y_phi phi_y2_phi) (b_0 phi_z_phi phi_z2_phi)) in - let m_beta = + let m_beta = let beta g = Mat.init_cols n_mo n_mo (fun i j -> (g.{i,i} -. g.{j,j})**2. -. 4. *. g.{i,j}**2.) in Mat.add (beta phi_x_phi) (Mat.add (beta phi_y_phi) (beta phi_z_phi)) - in - let m_gamma = + in + let m_gamma = let gamma g = Mat.init_cols n_mo n_mo (fun i j -> 4. *. g.{i,j} *. (g.{i,i} -. g.{j,j}) ) in Mat.add (gamma phi_x_phi) (Mat.add (gamma phi_y_phi) (gamma phi_z_phi)) - in + in let m_theta = Mat.init_cols n_mo n_mo (fun i j -> - if i = j + if i = j then 0. else +. 0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j}) 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 *) @@ -327,27 +327,27 @@ let localize mo_basis = let n_mo = Mat.dim2 m_C in let phi_x_phi = - Multipole.matrix_x multipoles + Multipole.matrix_x multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C - in + in let phi_y_phi = - Multipole.matrix_y multipoles + Multipole.matrix_y multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in let phi_z_phi = - Multipole.matrix_z multipoles + Multipole.matrix_z multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C - in + in - let m_b12= + let m_b12= let b12 g = Mat.init_cols n_mo 2 ( fun i j -> - if j = 1 + if j = 1 then (g.{i,i} -. g.{j1,j1}) *. g.{i,j1} else (g.{i,i} -. g.{j2,j2}) *. g.{i,j2}) - in + in Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi)) - in + in let m_a12 = let a12 g = Mat.init_cols n_mo 2 ( fun i j -> if j = 1 @@ -356,86 +356,82 @@ let localize mo_basis = in Mat.add (a12 phi_x_phi) ( Mat.add (a12 phi_y_phi) (a12 phi_z_phi)) in - (Mat.init_cols n_mo 2 ( fun i j -> - if i=j1 && j=1 + (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 if +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} >= 0. - then pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} - else -. pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j}), - 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.))) + 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 -> + (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)) 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 f2_theta_boys m_C j1 j2 = let multipoles = AOBasis.multipole ao_basis in let n_mo = Mat.dim2 m_C in let phi_x_phi = - Multipole.matrix_x multipoles + Multipole.matrix_x multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C - in + in let phi_y_phi = - Multipole.matrix_y multipoles + Multipole.matrix_y multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in let phi_z_phi = - Multipole.matrix_z multipoles + Multipole.matrix_z multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in - let phi_x2_phi = - Multipole.matrix_x2 multipoles - |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C - in - let phi_y2_phi = - Multipole.matrix_y2 multipoles + let phi_r2_phi = + Multipole.matrix_r2 multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in - let phi_z2_phi = - Multipole.matrix_z2 multipoles - |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C - in - - let m_beta = + + let m_beta = let beta g = Mat.init_cols n_mo 2 (fun i j -> if j = 1 then (g.{i,i} -. g.{j1,j1})**2. -. 4. *. g.{i,j1}**2. else (g.{i,i} -. g.{j2,j2})**2. -. 4. *. g.{i,j2}**2.) in Mat.add (beta phi_x_phi) (Mat.add (beta phi_y_phi) (beta phi_z_phi)) - in + in - let m_gamma = + 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 phi_x_phi) (Mat.add (gamma phi_y_phi) (gamma phi_z_phi)) - in + in let m_theta = Mat.init_cols n_mo 2 (fun i j -> - if i=j1 && j=1 + 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 + in let m_critere_B = Vec.init n_mo (fun i -> - phi_x2_phi.{i,i} -. phi_x_phi.{i,i}**2. +. phi_y2_phi.{i,i} -. phi_y_phi.{i,i}**2. +. phi_z2_phi.{i,i} -. phi_z_phi.{i,i}**2.) - + phi_r2_phi.{i,i} -. phi_x_phi.{i,i}**2. -. phi_y_phi.{i,i}**2. -. phi_z_phi.{i,i}**2.) + 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 *) + (* 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 m_alpha_d methode m_C = - let alpha methode = - match methode with + 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} @@ -446,12 +442,12 @@ let localize mo_basis = | "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 + 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 *) + (* 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 = @@ -468,13 +464,13 @@ let localize mo_basis = in (* 6. Norm of a matrix *) (* Matrix -> float *) - let norm m = - + let norm m = + let vec_m = Mat.as_vec m in let vec2 = Vec.sqr vec_m in sqrt(Vec.sum vec2) - in + 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 *) @@ -483,44 +479,44 @@ let localize mo_basis = 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.) + 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 - + in + (* Location of the max angle in the matrix *) (* angle matrix -> location of the max element *) (* matrix -> integer *) - let max_element alpha_m = + 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 + 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 alpha alpha_m = let i = indice_ii in let j = indice_jj in alpha_m.{i,j} in - let alpha_max = alpha alpha_m - + let alpha_max = alpha alpha_m + in if alpha_max < pi /. 2. then {alpha_max; indice_ii; indice_jj} @@ -532,11 +528,11 @@ let localize mo_basis = (* float -> matrix *) let f_R alpha = - Mat.init_cols 2 2 (fun i j -> - if i=j + Mat.init_cols 2 2 (fun i j -> + if i=j then cos alpha - else if i>j - then sin alpha + else if i>j + then sin alpha else -. sin alpha ) in @@ -544,7 +540,7 @@ let localize mo_basis = (* 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 *) @@ -553,15 +549,15 @@ let localize mo_basis = let m_angle angle= Mat.init_cols n_mo n_mo (fun i j -> - if i = j + if i = j then 0. - else if i>j + else if i>j then -. angle else angle) in - + let m = m_angle angle in - + (* Square of the matrix X² *) (* -> X . X *) (* -> matrix*) @@ -571,27 +567,27 @@ let localize mo_basis = (* 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 - + let tau = Vec.sqrt square_tau in + (* Calculation of cos tau from the vector tau *) (* -> cos(tau) *) (* -> integer *) 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} @@ -599,17 +595,17 @@ let localize mo_basis = in (* Calcul of sin(tau) from the vector tau *) - let sin_tau = Vec.sin tau in - + 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.) + if i=j + then sin_tau.{i} + else 0.) in (* Calculation of the transposed matrix of W (X²) *) @@ -630,12 +626,12 @@ let localize mo_basis = (* Calculation of the matrix tau^⁻1 *) (* -> matrix tau^-1 *) (* -> matrix *) - let m_tau_1 = + let m_tau_1 = - Mat.init_cols n_mo n_mo (fun i j -> - if i=j - then tau_1.{i} - else 0.) + Mat.init_cols n_mo n_mo (fun i j -> + if i=j + then tau_1.{i} + else 0.) in (* gemm mm (gemm m_cos_tau transpose_mm) -> W cos(tau) Wt *) @@ -649,7 +645,7 @@ let localize mo_basis = (* Sum of a + b -> Rotation matrix *) (* -> Rotation matrix *) (* -> matrix *) - let m_r = Mat.add a b + let m_r = Mat.add a b in gemm m_C m_r in @@ -658,19 +654,19 @@ let localize mo_basis = (* 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 + 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} + 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 + 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 *) @@ -679,12 +675,12 @@ let localize mo_basis = 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 + Mat.init_cols n_ao n_mo (fun i j -> + if j=indice_i then mat.{i,1} - else if j=indice_j + else if j=indice_j then mat.{i,2} else 0.) in @@ -693,14 +689,14 @@ let localize mo_basis = 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 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 + Mat.init_cols n_mo n_mo (fun i j -> + if j=indice_i then mat.{i,1} - else if j=indice_j + else if j=indice_j then mat.{i,2} else if i = indice_i then -. mat.{j,1} @@ -710,9 +706,9 @@ let localize mo_basis = 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 *) + (* 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 = + let f_interm m_C mat = Mat.sub m_C mat in @@ -722,7 +718,7 @@ let localize mo_basis = 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 n_rec_alpha = 10 in let alphaij = new_m_alpha m_alpha m_C n_rec_alpha in (* Value of the angle *) @@ -759,8 +755,8 @@ let localize mo_basis = (* 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 *) + + (* 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 *) @@ -768,7 +764,7 @@ let localize mo_basis = (* 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 @@ -792,27 +788,27 @@ let localize mo_basis = (* 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 *) + (* 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 + 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 @@ -823,8 +819,8 @@ let localize mo_basis = 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. + 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 @@ -833,7 +829,7 @@ let localize mo_basis = (* Function to create a integer list from a vector of float *) (* float vector -> integer list *) - let int_list vec = + let int_list vec = let float_list = Vec.to_list vec in let g a = int_of_float(a) @@ -844,7 +840,7 @@ let localize mo_basis = (* 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 list = + let miss_elem mat list = let n_mo = Mat.dim2 mat in let vec = Vec.init (n_mo) (fun i -> @@ -852,7 +848,7 @@ let localize mo_basis = then 0. else float_of_int(i)) in - let list_int = int_list vec + let list_int = int_list vec in List.filter (fun x -> x > 0) list_int in @@ -874,18 +870,18 @@ let localize mo_basis = (* Function to create a matrix from a list*) (* Coefficient matrix n_ao x n_mo, list of orbitals -> matrix with these orbitals *) - (* matrix, integer list -> matrix*) - let create_mat mat list = + (* matrix, integer list -> matrix*) + let create_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 = List.map f list + let vec_list = List.map f list in Mat.of_col_vecs_list vec_list in (* List of the occupied MOs *) - let list_occ = + let list_occ = let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i)) in @@ -893,11 +889,11 @@ let localize mo_basis = in (* List of the core MOs *) - let list_core = + let list_core = let vec_core = Vec.init (n_core) (fun i -> float_of_int(i)) in int_list vec_core in - + (* Function to create a list of pi MOs *) (* Coefficient matrix n_ao x n_mo -> list of pi orbitals *) (* matrix -> integer list *) @@ -906,7 +902,7 @@ let localize mo_basis = let m_pi = Mat.init_cols n_ao (n_mo+1) ( fun i j -> if j=1 && m_C.{i,j}**2. < 10e-8**2. then 0. - else if j=1 + else if j=1 then m_C.{i,j} else if mat.{i,j-1}**2. < 10e-8**2. then 0. @@ -928,7 +924,7 @@ let localize mo_basis = else var :: (remove x tail) in remove 0 list_pi - in + in (* Function to create a list of the (n-x_end) first MOs of a matrix *) (* Coefficient matrix n_ao x n_mo -> list of the (n-x_end) first MOs of the matrix *) @@ -949,15 +945,15 @@ let localize mo_basis = in int_list vec_12 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*) + (* 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 + let m = Array.append v_mat_1 v_mat_2 in Mat.of_col_vecs m in @@ -979,10 +975,10 @@ let localize mo_basis = (* 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 - (* Localization function *) + (* Localization function *) (* Coefficient matrix n_ao x n_mo -> localized coefficient matrix n_ao x n_mo *) (* matrix -> matrix *) let f_localise mat = @@ -996,40 +992,40 @@ let localize mo_basis = else rotate mat pre_angle in if Mat.dim2 mat = 0 - then new_mat + then new_mat else localise localisation new_mat methode epsilon iteration cc in - + (* Core MOs localization *) - let loc_m_core = - Printf.printf "%s\n" "I/III"; + let loc_m_core = + Printf.printf "%s\n" "I/III"; Printf.printf "%s\n" "Localization of core orbitals"; Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 m_core); Printf.printf "%s\n" ""; - + f_localise m_core in - + Printf.printf "%s\n" ""; Printf.printf "%s\n" "End I/III"; Printf.printf "%s\n" ""; (* Localization function without separation *) let f_loc_assemble_1 mat = - + Printf.printf "%s\n" ""; Printf.printf "%s\n" "Localization 1/1"; Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat); - + f_localise mat in - + Printf.printf "%s\n" ""; (* Localization function with 1 MOs separation *) let f_loc_assemble_12 mat list_12 = - - if List.length(list_12 mat) + List.length(miss_elem mat (list_12 mat)) <> (Mat.dim2 mat) + + if List.length(list_12 mat) + List.length(miss_elem mat (list_12 mat)) <> (Mat.dim2 mat) then Printf.printf "%s\n" "Bug : List length problem (list_12 or miss_elem list_12)"; - + let mat_1, mat_2 = split_mat mat (list_12 mat) in Printf.printf "%s\n" ""; @@ -1056,9 +1052,9 @@ let localize mo_basis = if List.length(list_2ab mat_2) + List.length(miss_elem mat_2 (list_2ab mat_2)) <> (Mat.dim2 mat_2) then Printf.printf "%s\n" "Bug : List length problem (list_2ab or miss_elem list_2ab)"; - + let mat_2a, mat_2b = split_mat mat_2 (list_2ab mat_2) in - + Printf.printf "%s\n" ""; Printf.printf "%s\n" "Localization 1/3"; Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_1); @@ -1077,7 +1073,7 @@ let localize mo_basis = let m_loc_2ab = assemble m_loc_2a m_loc_2b in - + assemble m_loc_1 m_loc_2ab in @@ -1128,7 +1124,7 @@ let localize mo_basis = (* Localization function with 4 separations "by hands" *) let f_loc_assemble_main mat list_1 list_2 list_3 list_4 = - + if (List.length(list_1) + List.length(list_2) + List.length(list_3) + List.length(list_4)) <> (Mat.dim2 mat) then Printf.printf "%s\n" "Bug : List length problem"; @@ -1136,7 +1132,7 @@ let localize mo_basis = let mat_2 = create_mat mat list_2 in let mat_3 = create_mat mat list_3 in let mat_4 = create_mat mat list_4 in - + Printf.printf "%s\n" ""; Printf.printf "%s\n" "Localization 1/4"; Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_1); @@ -1155,13 +1151,13 @@ let localize mo_basis = Printf.printf "%s\n" "Localization 4/4"; Printf.printf "%s %i\n" "Number of orbitals" (Mat.dim2 mat_4); - let loc_mat_4 = f_localise mat_4 - + let loc_mat_4 = f_localise mat_4 + in Printf.printf "%s\n" ""; assemble (assemble loc_mat_1 loc_mat_2) (assemble loc_mat_3 loc_mat_4) in - + (* Pattern matching for the separation type of the occupied MOs *) let m_assemble_occ, m_assemble_occu = @@ -1174,7 +1170,7 @@ let localize mo_basis = | "pi/sigma" -> let m_assemble_occu = f_loc_assemble_12 m_occu list_pi in m_assemble_occu | "pi/sigma/end" -> let m_assemble_occu = f_loc_assemble_123 m_occu list_pi list_x_end in m_assemble_occu | "pi/1/2" -> let m_assemble_occu = f_loc_assemble_123 m_occu list_pi list_s_12 in m_assemble_occu - | "1/2" -> let m_assemble_occu = f_loc_assemble_12 m_occu list_s_12 in m_assemble_occu + | "1/2" -> let m_assemble_occu = f_loc_assemble_12 m_occu list_s_12 in m_assemble_occu | "pi1/pi2/1/2" -> let m_assemble_occu = f_loc_assemble_1234 m_occu list_pi list_s_12 list_s_12 in m_assemble_occu | "By_hands" -> let m_assemble_occu = f_loc_assemble_main m_occu (miss_elem m_occ []) [] [] [] in m_assemble_occu | _ -> invalid_arg "Unknown separation type of valence orbitals" @@ -1186,7 +1182,7 @@ let localize mo_basis = (* Pattern matching for the separation type of the virtual MOs *) let m_assemble_vir = - + let m_assemble_vir = Printf.printf "%s\n" "III/III"; Printf.printf "%s\n" "Localization of virtual orbitals"; @@ -1215,41 +1211,41 @@ let localize mo_basis = (* 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= + let distrib mat= if Mat.dim1 mat = 0 - then Mat.as_vec mat + then Mat.as_vec mat else (let multipoles = AOBasis.multipole ao_basis in let n_mo = Mat.dim2 mat in let phi_x_phi = - Multipole.matrix_x multipoles + Multipole.matrix_x multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat - in + in let phi_y_phi = - Multipole.matrix_y multipoles + Multipole.matrix_y multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat in let phi_z_phi = - Multipole.matrix_z multipoles + Multipole.matrix_z multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat in let phi_x2_phi = - Multipole.matrix_x2 multipoles + Multipole.matrix_x2 multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat - in + in let phi_y2_phi = - Multipole.matrix_y2 multipoles + Multipole.matrix_y2 multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat in let phi_z2_phi = - Multipole.matrix_z2 multipoles + Multipole.matrix_z2 multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:mat - in + in Vec.init n_mo (fun i -> phi_x2_phi.{i,i} -. phi_x_phi.{i,i}**2. +. phi_y2_phi.{i,i} -. phi_y_phi.{i,i}**2. +. phi_z2_phi.{i,i} -. phi_z_phi.{i,i}**2.)) in - - (* Sorting function *) + + (* Sorting function *) (* tri insertion (<) list -> list in ascending order *) (* parameter, list -> list *) let rec insere comparaison elem liste = match liste with @@ -1266,13 +1262,13 @@ let localize mo_basis = (* 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 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 + if j=1 then vec_distrib.{i} else vec_tri_distrib.{i}) @@ -1295,14 +1291,14 @@ let localize mo_basis = in (* Function to compute the standard deviation *) - let f_stand_dev mat = + let f_stand_dev mat = if (Vec.dim (distrib mat)) = 0 then 0. - else sqrt(abs_float(f_variance mat)) + else sqrt(abs_float(f_variance mat)) in - + (* Fonction de calcul de la mediane *) - let f_median mat = + let f_median mat = if (Vec.dim (distrib mat)) = 0 then 0. else @@ -1311,13 +1307,13 @@ let localize mo_basis = 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. + 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 analyse mat = let average = f_average mat in let variance = f_variance mat in let stand_dev = f_stand_dev mat in @@ -1346,7 +1342,7 @@ let localize mo_basis = 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); @@ -1371,6 +1367,6 @@ let localize mo_basis = Printf.printf "%s %s\n" "Virtual orbitals after localization : " (analyse m_assemble_vir); Printf.printf "%s\n" ""; -m_assemble_loc +m_assemble_loc