From c4729da2a3215d2dd2bf079a80882dc6393f0866 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 10 Jul 2020 01:07:38 +0200 Subject: [PATCH] Optimization and cleaning --- MOBasis/Localisation.ml | 178 ++++++++++++++++++---------------------- 1 file changed, 81 insertions(+), 97 deletions(-) diff --git a/MOBasis/Localisation.ml b/MOBasis/Localisation.ml index c74449d..9985fd1 100644 --- a/MOBasis/Localisation.ml +++ b/MOBasis/Localisation.ml @@ -225,38 +225,48 @@ let localize mo_basis = let multipoles = AOBasis.multipole ao_basis in let n_mo = Mat.dim2 m_C in - let phi_x_phi = + let x = Multipole.matrix_x multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in - let phi_y_phi = + let y = Multipole.matrix_y multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in - let phi_z_phi = + let z = 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)) + 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 = - 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)) - + 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 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}) - ,0.5 *. 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 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 *) @@ -400,27 +410,24 @@ let localize mo_basis = let multipoles = AOBasis.multipole ao_basis in let n_mo = Mat.dim2 m_C in - let phi_x_phi = + let x = Multipole.matrix_x multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in - let phi_y_phi = + let y = Multipole.matrix_y multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in - let phi_z_phi = + let z = Multipole.matrix_z multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in - let phi_r2_phi = + let r2 = Multipole.matrix_r2 multipoles |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C in let m_beta = - let x = phi_x_phi - and y = phi_y_phi - and z = phi_z_phi in Mat.init_cols n_mo 2 (fun i j -> if j = 1 then @@ -448,7 +455,7 @@ let localize mo_basis = 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 Mat.add (gamma x) (Mat.add (gamma y) (gamma z)) in let m_theta = Mat.init_cols n_mo 2 (fun i j -> @@ -459,7 +466,7 @@ let localize mo_basis = else +. 0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j}) in let m_critere_B = Vec.init n_mo (fun i -> - phi_r2_phi.{i,i} -. phi_x_phi.{i,i}**2. -. phi_y_phi.{i,i}**2. -. phi_z_phi.{i,i}**2.) + r2.{i,i} -. x.{i,i}**2. -. y.{i,i}**2. -. z.{i,i}**2.) in m_theta, Vec.sum( m_critere_B) in @@ -517,17 +524,19 @@ let localize mo_basis = (* 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} ) + 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 *) @@ -558,9 +567,10 @@ let localize mo_basis = 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) + 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 *) @@ -572,7 +582,7 @@ let localize mo_basis = if i=j then cos alpha else if i>j - then sin alpha + then sin alpha else -. sin alpha ) in @@ -620,7 +630,7 @@ let localize mo_basis = (* Calculation of cos tau from the vector tau *) (* -> cos(tau) *) - (* -> integer *) + (* -> float *) let cos_tau = Vec.cos tau in (* Matrix cos tau *) @@ -648,11 +658,6 @@ let localize mo_basis = else 0.) in - (* Calculation of the transposed matrix of W (X²) *) - (* -> transposed matrix of W *) - (* -> matrix *) - let transpose_mm = Mat.transpose_copy mm in - (* Calculation of the vector tau^-1 *) (* -> vector tau^-1 *) (* -> vector *) @@ -674,17 +679,17 @@ let localize mo_basis = else 0.) in - (* gemm mm (gemm m_cos_tau transpose_mm) -> W cos(tau) Wt *) + (* W cos(tau) Wt *) (* -> matrix *) - let a = gemm mm (gemm m_cos_tau transpose_mm) in + let a = gemm mm @@ gemm ~transb:`T m_cos_tau mm in - (* gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) -> W tau^-1 sin(tau) Wt X *) + (* W tau^-1 sin(tau) Wt X *) (* -> matrix *) - let b = gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) in + 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 *) + (* -> matrix *) let m_r = Mat.add a b in gemm m_C m_r @@ -705,7 +710,7 @@ let localize mo_basis = (* 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 *) + (* matrix, matrix -> matrix *) let f_Ksi_tilde m_R m_Ksi = gemm m_Ksi m_R in @@ -867,41 +872,30 @@ let localize mo_basis = (* IV : Functions to split and assemble matrix *) - (* Function to create a integer list from a vector of float *) - (* float vector -> integer list *) - 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 - in (* 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 lst = 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 + 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 list = + 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 list in - let list_2 = miss_elem mat list 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) @@ -911,27 +905,23 @@ 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 = + let create_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 = List.map f list + let vec_list = List.map f lst in Mat.of_col_vecs_list vec_list in (* List of the occupied MOs *) let list_occ = - - let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i)) - in - int_list vec_occ + Util.list_range 1 nocc in (* List of the core MOs *) let list_core = - let vec_core = Vec.init (n_core) (fun i -> float_of_int(i)) - in int_list vec_core + Util.list_range 1 n_core in (* Function to create a list of pi MOs *) @@ -940,26 +930,24 @@ let localize mo_basis = 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 && m_C.{i,j}**2. < 10e-8**2. + if j=1 && abs_float m_C.{i,j} < 1.e-8 then 0. else if j=1 then m_C.{i,j} - else if mat.{i,j-1}**2. < 10e-8**2. + else if abs_float mat.{i,j-1} < 1.e-8 then 0. else mat.{i,j-1}) in let vec = Mat.to_col_vecs m_pi in - let vec_dot = Vec.init ((Mat.dim2 mat)) (fun i -> dot vec.(0) vec.(i)) in - let vec_pi = Vec.init ((Mat.dim2 mat)) (fun i -> - if vec_dot.{i} = 0. - then float_of_int(i) - else 0.) + (* 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 vec_dot.{i} = 0. then i else 0) in - let list_pi = int_list vec_pi in - let rec remove x list = - match list with - | [] -> [] - | var :: tail -> if var = x + let rec remove x = function + | [] -> [] + | var :: tail -> if var = x then remove x tail else var :: (remove x tail) in @@ -971,9 +959,7 @@ let localize mo_basis = (* matrix -> integer list *) let list_x_end mat = let n = Mat.dim2 mat in - let vec_x_end = Vec.init (n-x_end) (fun i -> float_of_int(i)) - in - int_list vec_x_end + Util.list_range 1 (n-x_end) in (* Function to create a list of the (n/2) first MOs of a matrix *) @@ -981,9 +967,7 @@ let localize mo_basis = (* matrix -> integer list *) let list_s_12 mat = let n = Mat.dim2 mat in - let vec_12 = Vec.init (n/2) (fun i -> float_of_int(i)) - in - int_list vec_12 + Util.list_range 1 (n/2) in (* Function to assemble to matrix with the same number of line *)