LCPQ
/
QCaml
mirror of https://gitlab.com/scemama/QCaml.git
10
1
Fork 0
Code Releases Activity

Optimization and cleaning

This commit is contained in:
Anthony Scemama 2020-07-10 01:07:38 +02:00
parent 5456852b0e
commit c4729da2a3
1 changed files with 81 additions and 97 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

178
MOBasis/Localisation.ml
View File

@ -225,38 +225,48 @@ let localize mo_basis =
let multipoles = AOBasis.multipole ao_basis in let multipoles = AOBasis.multipole ao_basis in
let n_mo = Mat.dim2 m_C in let n_mo = Mat.dim2 m_C in
let phi_x_phi = let x =
Multipole.matrix_x multipoles Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in in
let phi_y_phi = let y =
Multipole.matrix_y multipoles Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in in
let phi_z_phi = let z =
Multipole.matrix_z multipoles Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C |> 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 n_mo ( fun i j -> Mat.init_cols n_mo n_mo ( fun i j ->
(g.{i,i} -. g.{j,j}) *. g.{i,j}) (x.{i,i} -. x.{j,j}) *. x.{i,j} +.
in (y.{i,i} -. y.{j,j}) *. y.{i,j} +.
Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi)) (z.{i,i} -. z.{j,j}) *. z.{i,j})
in in
let m_a12 = let m_a12 =
let a12 g = Mat.init_cols n_mo n_mo ( fun i j -> 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}))) let x' = x.{i,i} -. x.{j,j}
in and y' = y.{i,i} -. y.{j,j}
Mat.add (a12 phi_x_phi) ( Mat.add (a12 phi_y_phi) (a12 phi_z_phi)) 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 in
(Mat.init_cols n_mo n_mo ( fun i j -> (Mat.init_cols n_mo n_mo ( fun i j ->
if i=j if i=j
then 0. then 0.
else if +. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} >= 0. else
then pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j} let w = atan2 m_b12.{i,j} m_a12.{i,j} in
else -. pi /. 4. -. 0.25 *. atan2 m_b12.{i,j} m_a12.{i,j}) if w >= 0.
,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.))) 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 in
(* 4. Function for the original Boys localization *) (* 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 *) (* 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 multipoles = AOBasis.multipole ao_basis in
let n_mo = Mat.dim2 m_C let n_mo = Mat.dim2 m_C
in in
let phi_x_phi = let x =
Multipole.matrix_x multipoles Multipole.matrix_x multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in in
let phi_y_phi = let y =
Multipole.matrix_y multipoles Multipole.matrix_y multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in in
let phi_z_phi = let z =
Multipole.matrix_z multipoles Multipole.matrix_z multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in in
let phi_r2_phi = let r2 =
Multipole.matrix_r2 multipoles Multipole.matrix_r2 multipoles
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
in in
let m_beta = 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 -> Mat.init_cols n_mo 2 (fun i j ->
if j = 1 then if j = 1 then
@ -448,7 +455,7 @@ let localize mo_basis =
if j = 1 if j = 1
then 4. *. g.{i,j1} *. (g.{i,i} -. g.{j1,j1}) then 4. *. g.{i,j1} *. (g.{i,i} -. g.{j1,j1})
else 4. *. g.{i,j2} *. (g.{i,i} -. g.{j2,j2})) 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 in
let m_theta = Mat.init_cols n_mo 2 (fun i j -> 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}) else +. 0.25 *. atan2 m_gamma.{i,j} m_beta.{i,j})
in in
let m_critere_B = Vec.init n_mo (fun i -> 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 m_theta, Vec.sum( m_critere_B)
in in
@ -517,17 +524,19 @@ let localize mo_basis =
(* matrix, matrix, integer -> float, integer, integer *) (* matrix, matrix, integer -> float, integer, integer *)
let rec new_m_alpha m_alpha m_C n_rec_alpha= 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 n_mo = Mat.dim2 m_C in
let alpha_m = let alpha_m =
if n_rec_alpha == 0 if n_rec_alpha == 0 then
then m_alpha m_alpha
else else
Mat.init_cols n_mo n_mo (fun i j -> Mat.init_cols n_mo n_mo (fun i j ->
if (m_alpha.{i,j}) >= (pi /. 2.) if m_alpha.{i,j} >= pi_2 then
then (m_alpha.{i,j} -. ( pi /. 2.)) m_alpha.{i,j} -. pi_2
else if m_alpha.{i,j} <= -. pi /. 2. else if m_alpha.{i,j} <= -. pi_2 then
then (m_alpha.{i,j} +. ( pi /. 2.)) m_alpha.{i,j} +. pi_2
else m_alpha.{i,j} ) else m_alpha.{i,j}
)
in in
(* Location of the max angle in the matrix *) (* Location of the max angle in the matrix *)
@ -558,9 +567,10 @@ let localize mo_basis =
let alpha_max = alpha alpha_m let alpha_max = alpha alpha_m
in in
if alpha_max < pi /. 2. if alpha_max < pi_2 then
then {alpha_max; indice_ii; indice_jj} {alpha_max; indice_ii; indice_jj}
else new_m_alpha alpha_m m_C (n_rec_alpha-1) else
new_m_alpha alpha_m m_C (n_rec_alpha-1)
in in
(* 8. Rotation matrix 2x2 *) (* 8. Rotation matrix 2x2 *)
@ -572,7 +582,7 @@ let localize mo_basis =
if i=j if i=j
then cos alpha then cos alpha
else if i>j else if i>j
then sin alpha then sin alpha
else -. sin alpha ) else -. sin alpha )
in in
@ -620,7 +630,7 @@ let localize mo_basis =
(* Calculation of cos tau from the vector tau *) (* Calculation of cos tau from the vector tau *)
(* -> cos(tau) *) (* -> cos(tau) *)
(* -> integer *) (* -> float *)
let cos_tau = Vec.cos tau in let cos_tau = Vec.cos tau in
(* Matrix cos tau *) (* Matrix cos tau *)
@ -648,11 +658,6 @@ let localize mo_basis =
else 0.) else 0.)
in 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 *) (* Calculation of the vector tau^-1 *)
(* -> vector tau^-1 *) (* -> vector tau^-1 *)
(* -> vector *) (* -> vector *)
@ -674,17 +679,17 @@ let localize mo_basis =
else 0.) else 0.)
in in
(* gemm mm (gemm m_cos_tau transpose_mm) -> W cos(tau) Wt *) (* W cos(tau) Wt *)
(* -> matrix *) (* -> 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 *) (* -> 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 *) (* Sum of a + b -> Rotation matrix *)
(* -> Rotation matrix *) (* -> Rotation matrix *)
(* -> matrix *) (* -> matrix *)
let m_r = Mat.add a b let m_r = Mat.add a b
in in
gemm m_C m_r 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 *) (* 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 *) (* 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 let f_Ksi_tilde m_R m_Ksi = gemm m_Ksi m_R
in in
@ -867,41 +872,30 @@ let localize mo_basis =
(* IV : Functions to split and assemble matrix *) (* 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*) (* 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 *) (* Coefficient matrix n_ao x n_mo, list -> list of missing element of the previous list *)
(* matrix, integer list -> integer list *) (* matrix, integer list -> integer list *)
let miss_elem mat list = let miss_elem mat lst =
let n_mo = Mat.dim2 mat in let n_mo = Mat.dim2 mat in
let vec = Vec.init (n_mo) (fun i -> Util.list_range 1 n_mo
if List.mem i list |> List.map (fun i ->
then 0. if List.mem i lst
else float_of_int(i)) then 0
in else i)
let list_int = int_list vec |> List.filter (fun x -> x > 0)
in
List.filter (fun x -> x > 0) list_int
in 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 *) (* 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 *) (* Coefficient matrix n_ao x n_mo, list -> matrix, matrix *)
(* matrix, integer 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 vec_of_mat = Mat.to_col_vecs mat in
let f a = vec_of_mat.(a-1) in let f a = vec_of_mat.(a-1) in
let vec_list_1 = List.map f list in let vec_list_1 = List.map f lst in
let list_2 = miss_elem mat list in let list_2 = miss_elem mat lst in
let vec_list_2 = List.map f list_2 let vec_list_2 = List.map f list_2
in in
(Mat.of_col_vecs_list vec_list_1,Mat.of_col_vecs_list vec_list_2) (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*) (* Function to create a matrix from a list*)
(* Coefficient matrix n_ao x n_mo, list of orbitals -> matrix with these orbitals *) (* Coefficient matrix n_ao x n_mo, list of orbitals -> matrix with these orbitals *)
(* matrix, integer list -> matrix*) (* 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 vec_of_mat = Mat.to_col_vecs mat in
let f a = vec_of_mat.(a-1) 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 in
Mat.of_col_vecs_list vec_list Mat.of_col_vecs_list vec_list
in in
(* List of the occupied MOs *) (* List of the occupied MOs *)
let list_occ = let list_occ =
Util.list_range 1 nocc
let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))
in
int_list vec_occ
in in
(* List of the core MOs *) (* List of the core MOs *)
let list_core = let list_core =
let vec_core = Vec.init (n_core) (fun i -> float_of_int(i)) Util.list_range 1 n_core
in int_list vec_core
in in
(* Function to create a list of pi MOs *) (* Function to create a list of pi MOs *)
@ -940,26 +930,24 @@ let localize mo_basis =
let list_pi mat = let list_pi mat =
let n_mo = Mat.dim2 mat in let n_mo = Mat.dim2 mat in
let m_pi = Mat.init_cols n_ao (n_mo+1) ( fun i j -> 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. then 0.
else if j=1 else if j=1
then m_C.{i,j} 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. then 0.
else mat.{i,j-1}) else mat.{i,j-1})
in in
let vec = Mat.to_col_vecs m_pi 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 (* TODO gemv *)
let vec_pi = Vec.init ((Mat.dim2 mat)) (fun i -> let vec_dot = Vec.init (Mat.dim2 mat) (fun i -> dot vec.(0) vec.(i)) in
if vec_dot.{i} = 0. let list_pi =
then float_of_int(i) Util.list_range 1 (Mat.dim2 mat)
else 0.) |> List.map (fun i -> if vec_dot.{i} = 0. then i else 0)
in in
let list_pi = int_list vec_pi in let rec remove x = function
let rec remove x list = | [] -> []
match list with | var :: tail -> if var = x
| [] -> []
| var :: tail -> if var = x
then remove x tail then remove x tail
else var :: (remove x tail) else var :: (remove x tail)
in in
@ -971,9 +959,7 @@ let localize mo_basis =
(* matrix -> integer list *) (* matrix -> integer list *)
let list_x_end mat = let list_x_end mat =
let n = Mat.dim2 mat in let n = Mat.dim2 mat in
let vec_x_end = Vec.init (n-x_end) (fun i -> float_of_int(i)) Util.list_range 1 (n-x_end)
in
int_list vec_x_end
in in
(* Function to create a list of the (n/2) first MOs of a matrix *) (* 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 *) (* matrix -> integer list *)
let list_s_12 mat = let list_s_12 mat =
let n = Mat.dim2 mat in let n = Mat.dim2 mat in
let vec_12 = Vec.init (n/2) (fun i -> float_of_int(i)) Util.list_range 1 (n/2)
in
int_list vec_12
in in
(* Function to assemble to matrix with the same number of line *) (* Function to assemble to matrix with the same number of line *)