908 lines
25 KiB
OCaml
908 lines
25 KiB
OCaml
let png_image = print_endline ;;
|
|
|
|
(* --------- *)
|
|
|
|
|
|
|
|
|
|
|
|
(* --------- *)
|
|
|
|
open Lacaml.D
|
|
|
|
|
|
|
|
let charge = 0
|
|
|
|
let xyz_string =
|
|
"3
|
|
Water
|
|
O 0. 0. 0.
|
|
H -0.756950272703377558 0. -0.585882234512562827
|
|
H 0.756950272703377558 0. -0.585882234512562827
|
|
"
|
|
|
|
(* Fonction création chaîne linéaire de n H *)
|
|
|
|
let xyz d n =
|
|
let accu = ""
|
|
in
|
|
let rec toto accu d n =
|
|
let accu =
|
|
if n=0
|
|
then accu ^ ""
|
|
else
|
|
match n with
|
|
| 1 -> "H 0. 0. 0.\n" ^ accu
|
|
| x -> toto ("H" ^ " " ^ string_of_float( d *. float_of_int(n-1)) ^ " 0. 0.\n" ^ accu ) d (n-1)
|
|
in
|
|
accu
|
|
in string_of_int(n) ^ "\nH" ^ string_of_int(n) ^ "\n" ^ toto accu d n;;
|
|
|
|
let xyz_string = xyz 1.8 6;;
|
|
|
|
let charge=1;;
|
|
|
|
let xyz_string =
|
|
"28
|
|
Cyanine C11
|
|
C 0.000000 0.000000 0.553396
|
|
H 0.000000 0.000000 1.639052
|
|
C 0.000000 1.212451 -0.115894
|
|
C 0.000000 -1.212451 -0.115894
|
|
H 0.000000 1.184500 -1.203671
|
|
H 0.000000 -1.184500 -1.203671
|
|
C 0.000000 2.463852 0.491562
|
|
C 0.000000 -2.463852 0.491562
|
|
H 0.000000 2.518814 1.575709
|
|
H 0.000000 -2.518814 1.575709
|
|
C 0.000000 3.632731 -0.239322
|
|
C 0.000000 -3.632731 -0.239322
|
|
H 0.000000 3.543995 -1.323893
|
|
H 0.000000 -3.543995 -1.323893
|
|
C 0.000000 4.922617 0.294909
|
|
C 0.000000 -4.922617 0.294909
|
|
H 0.000000 5.053656 1.372148
|
|
H 0.000000 -5.053656 1.372148
|
|
C 0.000000 6.023589 -0.522235
|
|
C 0.000000 -6.023589 -0.522235
|
|
H 0.000000 5.885105 -1.598471
|
|
H 0.000000 -5.885105 -1.598471
|
|
N 0.000000 7.289195 -0.119132
|
|
N 0.000000 -7.289195 -0.119132
|
|
H 0.000000 7.528514 0.857192
|
|
H 0.000000 -7.528514 0.857192
|
|
H 0.000000 8.044630 -0.778923
|
|
H 0.000000 -8.044630 -0.778923
|
|
"
|
|
|
|
let xyz_string =
|
|
"8
|
|
C1
|
|
C 0.000000 0.000000 0.424165
|
|
H 0.000000 0.000000 1.508704
|
|
N 0.000000 1.158087 -0.174215
|
|
N 0.000000 -1.158087 -0.174215
|
|
H 0.000000 1.258360 -1.178487
|
|
H 0.000000 2.005112 0.371143
|
|
H 0.000000 -2.005112 0.371143
|
|
H 0.000000 -1.258360 -1.178487
|
|
"
|
|
|
|
|
|
let basis_string =
|
|
"
|
|
HYDROGEN
|
|
S 4
|
|
1 1.301000E+01 1.968500E-02
|
|
2 1.962000E+00 1.379770E-01
|
|
3 4.446000E-01 4.781480E-01
|
|
4 1.220000E-01 5.012400E-01
|
|
S 1
|
|
1 1.220000E-01 1.000000E+00
|
|
P 1
|
|
1 7.270000E-01 1.0000000
|
|
|
|
CARBON
|
|
S 9
|
|
1 6.665000E+03 6.920000E-04
|
|
2 1.000000E+03 5.329000E-03
|
|
3 2.280000E+02 2.707700E-02
|
|
4 6.471000E+01 1.017180E-01
|
|
5 2.106000E+01 2.747400E-01
|
|
6 7.495000E+00 4.485640E-01
|
|
7 2.797000E+00 2.850740E-01
|
|
8 5.215000E-01 1.520400E-02
|
|
9 1.596000E-01 -3.191000E-03
|
|
S 9
|
|
1 6.665000E+03 -1.460000E-04
|
|
2 1.000000E+03 -1.154000E-03
|
|
3 2.280000E+02 -5.725000E-03
|
|
4 6.471000E+01 -2.331200E-02
|
|
5 2.106000E+01 -6.395500E-02
|
|
6 7.495000E+00 -1.499810E-01
|
|
7 2.797000E+00 -1.272620E-01
|
|
8 5.215000E-01 5.445290E-01
|
|
9 1.596000E-01 5.804960E-01
|
|
S 1
|
|
1 1.596000E-01 1.000000E+00
|
|
P 4
|
|
1 9.439000E+00 3.810900E-02
|
|
2 2.002000E+00 2.094800E-01
|
|
3 5.456000E-01 5.085570E-01
|
|
4 1.517000E-01 4.688420E-01
|
|
P 1
|
|
1 1.517000E-01 1.000000E+00
|
|
D 1
|
|
1 5.500000E-01 1.0000000
|
|
|
|
NITROGEN
|
|
S 9
|
|
1 9.046000E+03 7.000000E-04
|
|
2 1.357000E+03 5.389000E-03
|
|
3 3.093000E+02 2.740600E-02
|
|
4 8.773000E+01 1.032070E-01
|
|
5 2.856000E+01 2.787230E-01
|
|
6 1.021000E+01 4.485400E-01
|
|
7 3.838000E+00 2.782380E-01
|
|
8 7.466000E-01 1.544000E-02
|
|
9 2.248000E-01 -2.864000E-03
|
|
S 9
|
|
1 9.046000E+03 -1.530000E-04
|
|
2 1.357000E+03 -1.208000E-03
|
|
3 3.093000E+02 -5.992000E-03
|
|
4 8.773000E+01 -2.454400E-02
|
|
5 2.856000E+01 -6.745900E-02
|
|
6 1.021000E+01 -1.580780E-01
|
|
7 3.838000E+00 -1.218310E-01
|
|
8 7.466000E-01 5.490030E-01
|
|
9 2.248000E-01 5.788150E-01
|
|
S 1
|
|
1 2.248000E-01 1.000000E+00
|
|
P 4
|
|
1 1.355000E+01 3.991900E-02
|
|
2 2.917000E+00 2.171690E-01
|
|
3 7.973000E-01 5.103190E-01
|
|
4 2.185000E-01 4.622140E-01
|
|
P 1
|
|
1 2.185000E-01 1.000000E+00
|
|
D 1
|
|
1 8.170000E-01 1.0000000
|
|
|
|
"
|
|
|
|
let nuclei =
|
|
Nuclei.of_xyz_string xyz_string
|
|
|
|
let basis =
|
|
Basis.of_nuclei_and_basis_string nuclei basis_string
|
|
|
|
let simulation =
|
|
Simulation.make ~charge ~multiplicity:1 ~nuclei basis
|
|
|
|
let ao_basis =
|
|
Simulation.ao_basis simulation
|
|
|
|
let nocc =
|
|
let elec = Simulation.electrons simulation in
|
|
Electrons.n_alfa elec
|
|
|
|
|
|
|
|
|
|
let hf = HartreeFock.make simulation ;;
|
|
|
|
let ppf = Printing.ppf_dev_null ;;
|
|
Format.fprintf ppf "@[%a@]@." HartreeFock.pp hf ;;
|
|
|
|
let mo_basis = MOBasis.of_hartree_fock hf
|
|
|
|
let mo_coef = MOBasis.mo_coef mo_basis
|
|
|
|
|
|
(*
|
|
let m_X =
|
|
AOBasis.ortho ao_basis
|
|
*)
|
|
|
|
(*
|
|
let c_of_h m_H =
|
|
(* On exprime H dans la base orthonormale *)
|
|
let m_Hmo =
|
|
Util.xt_o_x m_H m_X (* H_mo = X^t H X *)
|
|
in
|
|
|
|
(* On diagonalise cet Hamiltonien *)
|
|
let m_C', _ =
|
|
Util.diagonalize_symm m_Hmo
|
|
in
|
|
|
|
(* On re-exprime les MOs dans la base des AOs (non-orthonormales) *)
|
|
gemm m_X m_C' (* C = X.C' *)
|
|
|
|
|
|
let m_C =
|
|
match Guess.make ~nocc ~guess:`Hcore ao_basis with
|
|
| Hcore m_H -> c_of_h m_H
|
|
| _ -> assert false
|
|
*)
|
|
|
|
(*
|
|
let m_P =
|
|
(* P = 2 C.C^t *)
|
|
gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C
|
|
*)
|
|
|
|
(*
|
|
let m_Hc, m_J, m_K =
|
|
let f =
|
|
Fock.make_rhf ~density:m_P ao_basis
|
|
in
|
|
Fock.(core f, coulomb f, exchange f)
|
|
*)
|
|
|
|
(*
|
|
let m_F =
|
|
Mat.add m_Hc (Mat.sub m_J m_K)
|
|
*)
|
|
|
|
(*let m_C =
|
|
let m_C', _ =
|
|
Util.xt_o_x m_F m_X
|
|
|> Util.diagonalize_symm
|
|
in
|
|
gemm m_X m_C'
|
|
*)
|
|
|
|
(*
|
|
let energy =
|
|
(Simulation.nuclear_repulsion simulation) +. 0.5 *.
|
|
Mat.gemm_trace m_P (Mat.add m_Hc m_F)
|
|
*)
|
|
|
|
(*
|
|
let rec iteration m_C n =
|
|
let m_P =
|
|
(* P = 2 C.C^t *)
|
|
gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C
|
|
in
|
|
|
|
let m_Hc, m_J, m_K =
|
|
let f =
|
|
Fock.make_rhf ~density:m_P ao_basis
|
|
in
|
|
Fock.(core f, coulomb f, exchange f)
|
|
in
|
|
|
|
let m_F =
|
|
Mat.add m_Hc (Mat.sub m_J m_K)
|
|
in
|
|
|
|
let m_C =
|
|
let m_C', _ =
|
|
Util.xt_o_x m_F m_X
|
|
|> Util.diagonalize_symm
|
|
in
|
|
gemm m_X m_C'
|
|
in
|
|
|
|
let energy =
|
|
(Simulation.nuclear_repulsion simulation) +. 0.5 *.
|
|
Mat.gemm_trace m_P (Mat.add m_Hc m_F)
|
|
in
|
|
Printf.printf "%f\n%!" energy;
|
|
if n > 0 then
|
|
iteration m_C (n-1)
|
|
|
|
let () = iteration m_C 20
|
|
*)
|
|
|
|
(* Construction de la matrice de rotation R de taille n par n *)
|
|
|
|
|
|
(* Définitions de base nécessaire pour la suite *)
|
|
let ee_ints = AOBasis.ee_ints ao_basis;;
|
|
let m_C = MOBasis.mo_coef mo_basis;;
|
|
let n_ao = Mat.dim1 m_C ;;
|
|
let n_mo = Mat.dim2 m_C ;;
|
|
let multipoles =
|
|
AOBasis.multipole ao_basis;;
|
|
|
|
let sum a =
|
|
Array.fold_left (fun accu x -> accu +. x) 0. a
|
|
|
|
let pi = 3.14159265358979323846264338;;
|
|
|
|
|
|
|
|
|
|
(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)
|
|
let f_alpha m_C =
|
|
|
|
let n_mo = Mat.dim2 m_C in
|
|
(*let t0 = Sys.time () in*)
|
|
|
|
let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in
|
|
let m_a12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in
|
|
let v_d = Vec.init n_mo (fun i -> 0.) in
|
|
|
|
(* Tableaux temporaires *)
|
|
let m_pqr =
|
|
Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)
|
|
in
|
|
let m_qr_i = Mat.create (n_ao*n_ao) n_mo in
|
|
let m_ri_j = Mat.create (n_ao*n_mo) n_mo in
|
|
let m_ij_k = Mat.create (n_mo*n_mo) n_mo in
|
|
|
|
Array.iter (fun s ->
|
|
(* Grosse boucle externe sur s *)
|
|
Array.iter (fun r ->
|
|
Array.iter (fun q ->
|
|
Array.iter (fun p ->
|
|
m_pqr.{p,q,r} <- ERI.get_phys ee_ints p q r s
|
|
) (Util.array_range 1 n_ao)
|
|
) (Util.array_range 1 n_ao)
|
|
) (Util.array_range 1 n_ao);
|
|
|
|
(* Conversion d'un tableau a 3 indices en une matrice nao x nao^2 *)
|
|
let m_p_qr =
|
|
Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]
|
|
|> Bigarray.array2_of_genarray
|
|
in
|
|
|
|
let m_qr_i =
|
|
(* (qr,i) = <i r|q s> = \sum_p <p r | q s> C_{pi} *)
|
|
gemm ~transa:`T ~c:m_qr_i m_p_qr m_C
|
|
in
|
|
|
|
let m_q_ri =
|
|
(* Transformation de la matrice (qr,i) en (q,ri) *)
|
|
Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo)
|
|
in
|
|
|
|
let m_ri_j =
|
|
(* (ri,j) = <i r | j s> = \sum_q <i r | q s> C_{bj} *)
|
|
gemm ~transa:`T ~c:m_ri_j m_q_ri m_C
|
|
in
|
|
|
|
let m_r_ij =
|
|
(* Transformation de la matrice (ri,j) en (r,ij) *)
|
|
Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo)
|
|
in
|
|
|
|
let m_ij_k =
|
|
(* (ij,k) = <i k | j s> = \sum_r <i r | j s> C_{rk} *)
|
|
gemm ~transa:`T ~c:m_ij_k m_r_ij m_C
|
|
in
|
|
|
|
let m_ijk =
|
|
(* Transformation de la matrice (ei,j) en (e,ij) *)
|
|
Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |]
|
|
|> Bigarray.array3_of_genarray
|
|
in
|
|
|
|
Array.iter (fun j ->
|
|
Array.iter (fun i ->
|
|
m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});
|
|
m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j} -.
|
|
0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_C.{s,i} +.
|
|
(m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_C.{s,j})
|
|
) (Util.array_range 1 n_mo);
|
|
v_d.{j} <- v_d.{j} +. m_ijk.{j,j,j} *. m_C.{s,j}
|
|
) (Util.array_range 1 n_mo)
|
|
) (Util.array_range 1 n_ao);
|
|
|
|
(*let t1 = Sys.time () in
|
|
Printf.printf "t = %f s\n%!" (t1 -. t0);*)
|
|
|
|
(Mat.init_cols n_mo n_mo ( fun i j ->
|
|
if i= j then 0.
|
|
else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2. ))))
|
|
),Vec.sum v_d);;
|
|
|
|
(*********************)
|
|
(*
|
|
f_alpha m_C;;
|
|
let m_alpha , s_D = f_alpha m_C;;
|
|
*)
|
|
|
|
(* Fonction calcul alpha Boys *)
|
|
let f_alpha_boys m_C =
|
|
let n_mo = Mat.dim2 m_C
|
|
in
|
|
let phi_x_phi =
|
|
Multipole.matrix_x multipoles
|
|
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
|
|
in
|
|
let phi_y_phi =
|
|
Multipole.matrix_y multipoles
|
|
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
|
|
in
|
|
let phi_z_phi =
|
|
Multipole.matrix_z multipoles
|
|
|> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C
|
|
in
|
|
|
|
let m_b12=
|
|
let b12 g = Mat.init_cols n_mo n_mo ( fun i j ->
|
|
(g.{i,i} -. g.{j,j}) *. g.{i,j})
|
|
|
|
in
|
|
Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi))
|
|
in
|
|
let m_a12 =
|
|
let a12 g = Mat.init_cols n_mo n_mo ( fun i j ->
|
|
g.{i,j} *. g.{i,j} -. 0.25 *. ((g.{i,i} -. g.{j,j}) *. (g.{i,i} -. g.{j,j})))
|
|
in
|
|
Mat.add (a12 phi_x_phi) ( Mat.add (a12 phi_y_phi) (a12 phi_z_phi))
|
|
in
|
|
(Mat.init_cols n_mo n_mo ( fun i j ->
|
|
if i=j
|
|
then 0.
|
|
else 0.25 *. acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) ))),
|
|
Vec.sum(Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)));;
|
|
|
|
(****************************)
|
|
(*
|
|
f_alpha_boys m_C;;
|
|
*)
|
|
|
|
(* Test méthode de calcul de alpha et de D ensemble *)
|
|
type alphad = {
|
|
m_alpha : Mat.t;
|
|
d : float;
|
|
}
|
|
|
|
let m_alpha_d methode m_C =
|
|
let alpha methode =
|
|
match methode with
|
|
| "Boys"
|
|
| "boys" -> let alpha_boys , d_boys = f_alpha_boys m_C in
|
|
{m_alpha = alpha_boys; d = d_boys}
|
|
| "ER"
|
|
| "er" -> let alpha_er , d_er = f_alpha m_C in
|
|
{m_alpha = alpha_er; d = d_er}
|
|
| _ -> invalid_arg "Unknown method, please enter Boys or ER"
|
|
|
|
in
|
|
alpha methode;;
|
|
|
|
(*************************)
|
|
(*
|
|
m_alpha_d "ER" ;;
|
|
m_alpha_d "Boys" ;;
|
|
let methode = "ER";;
|
|
let alphad = m_alpha_d methode m_C;;
|
|
let m_alpha = alphad.m_alpha;;
|
|
let d = alphad.d;;
|
|
*)
|
|
|
|
(* Test norme de alpha *)
|
|
|
|
let norme m =
|
|
let vec_m = Mat.as_vec m
|
|
in
|
|
let vec2 = Vec.sqr vec_m
|
|
in sqrt(Vec.sum vec2);;
|
|
|
|
(************************)
|
|
(*
|
|
norme_alpha m_alpha;;
|
|
*)
|
|
|
|
type alphaij = {
|
|
alpha_max : float;
|
|
indice_ii : int;
|
|
indice_jj : int;};;
|
|
|
|
(* Détermination alpha_max et ses indices i et j.
|
|
Si alpha max > pi/2 on soustrait pi/2 à la matrice des alphas de manière récursive *)
|
|
let rec new_m_alpha m_alpha m_C n_rec_alpha=
|
|
|
|
let n_mo = Mat.dim2 m_C
|
|
in
|
|
let alpha_m =
|
|
|
|
(*Printf.printf "%i\n%!" n_rec_alpha;*)
|
|
|
|
if n_rec_alpha == 0
|
|
then m_alpha
|
|
else Mat.init_cols n_mo n_mo (fun i j ->
|
|
if (m_alpha.{i,j}) > (pi /. 2.)
|
|
then (m_alpha.{i,j} -. ( pi /. 2.))
|
|
else if m_alpha.{i,j} < -. pi /. 2.
|
|
then (m_alpha.{i,j} +. ( pi /. 2.))
|
|
else if m_alpha.{i,j} < 0.
|
|
then -. m_alpha.{i,j}
|
|
else m_alpha.{i,j} )
|
|
in
|
|
|
|
(*Util.debug_matrix "alpha_m" alpha_m;*)
|
|
|
|
(* Détermination de l'emplacement du alpha max *)
|
|
let max_element3 alpha_m =
|
|
Mat.as_vec alpha_m
|
|
|> iamax
|
|
in
|
|
|
|
(* indice i et j du alpha max *)
|
|
let indice_ii, indice_jj =
|
|
let max = max_element3 alpha_m
|
|
in
|
|
(max - 1) mod n_mo +1, (max - 1) / n_mo +1
|
|
in
|
|
|
|
(* Valeur du alpha max*)
|
|
let alpha alpha_m =
|
|
let i = indice_ii
|
|
in
|
|
let j = indice_jj
|
|
in
|
|
|
|
(*Printf.printf "%i %i\n%!" i j;*)
|
|
|
|
alpha_m.{i,j}
|
|
|
|
in
|
|
let alpha_max = alpha alpha_m
|
|
in
|
|
|
|
(*Printf.printf "%f\n%!" alpha_max;*)
|
|
|
|
if alpha_max < pi /. 2.
|
|
then {alpha_max; indice_ii; indice_jj}
|
|
else new_m_alpha alpha_m m_C (n_rec_alpha-1);;
|
|
|
|
(*************************)
|
|
(*
|
|
let m_alpha,d = f_alpha m_C
|
|
let alphaij = new_m_alpha m_alpha m_C 3;;
|
|
alphaij.alpha_max;;
|
|
*)
|
|
|
|
(* Fonction de pattern matching pour localiser ou délocaliser *)
|
|
let alpha_v loc_deloc alphaij =
|
|
let alpha_loc = alphaij.alpha_max
|
|
in
|
|
let alpha_deloc = alphaij.alpha_max +. (pi /. 4.)
|
|
in
|
|
let choice loc_deloc =
|
|
match loc_deloc with
|
|
|"loc" -> alpha_loc
|
|
|"deloc" -> alpha_deloc
|
|
| _ -> invalid_arg "Unknown method, please enter loc or deloc"
|
|
in choice loc_deloc ;;
|
|
|
|
(* Matrice de rotation 2 par 2 *)
|
|
let f_R alpha =
|
|
Mat.init_cols 2 2 (fun i j ->
|
|
if i=j
|
|
then cos alpha
|
|
else if i>j
|
|
then sin alpha
|
|
else -. sin alpha )
|
|
;;
|
|
(*************************)
|
|
(*
|
|
let alpha = alphaij.alpha_max;; (* Fonction -> constante *)
|
|
f_R alpha;;
|
|
*)
|
|
(*
|
|
alpha_v "deloc" alphaij;;
|
|
let alpha = (alpha_v "loc" alphaij) ;;
|
|
f_R alpha ;;
|
|
*)
|
|
|
|
(*Uniquement pour pouvoir tester les fonctions après cette cellules*)
|
|
(*
|
|
(* Indice i et j du alpha max après calcul *)
|
|
let indice_i = alphaij.indice_ii;;
|
|
let indice_j = alphaij.indice_jj;;
|
|
let m_R = f_R alpha;;
|
|
*)
|
|
|
|
(* Fonction d'extraction des 2 vecteurs propres i et j de la matrice des OMs pour les mettres dans la matrice Ksi (n par 2)
|
|
pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)
|
|
let f_Ksi indice_i indice_j m_C =
|
|
let n_ao = Mat.dim1 m_C
|
|
in
|
|
Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_C.{i,indice_i} else m_C.{i,indice_j} );;
|
|
|
|
(*************************)
|
|
(*
|
|
let m_Ksi = f_Ksi indice_i indice_j m_C;;
|
|
*)
|
|
|
|
(* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle
|
|
on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)
|
|
let f_Ksi_tilde m_R m_Ksi m_C = gemm m_Ksi m_R;;
|
|
|
|
(*************************)
|
|
(*
|
|
let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi;;
|
|
*)
|
|
|
|
(* Fonction pour la création de matrice intermédiaire pour supprimer i j et ajouter i~ et j~ *)
|
|
let f_k mat indice_i indice_j m_C =
|
|
let n_mo = Mat.dim2 m_C
|
|
in
|
|
let n_ao = Mat.dim1 m_C
|
|
in
|
|
Mat.init_cols n_ao n_mo (fun i j ->
|
|
if j=indice_i
|
|
then mat.{i,1}
|
|
else if j=indice_j
|
|
then mat.{i,2}
|
|
else 0.)
|
|
(*************************)
|
|
(*
|
|
let m_Psi = f_Psi m_Ksi indice_i indice_j;;
|
|
let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j;;
|
|
*)
|
|
|
|
(* 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 *)
|
|
let f_interm m_C m_Psi = Mat.sub m_C m_Psi;;
|
|
|
|
(*************************)
|
|
(*
|
|
let m_interm = f_interm m_C m_Psi;;
|
|
let new_m_C m_C= Mat.add m_Psi_tilde m_interm;;
|
|
let m_new_m_C = new_m_C m_C;;
|
|
*)
|
|
|
|
(* Test localisation matrice rectangulaire et partielle *)
|
|
(*let toto = [4];;
|
|
let occ_m_C m_C toto= Mat.init_cols 4 3 (fun i j ->
|
|
if not (List.mem j toto)
|
|
then m_C.{i,j}
|
|
else 0.);;
|
|
|
|
let occ = occ_m_C m_C toto;;
|
|
*)
|
|
|
|
(* Localisation de Edminstion ou de Boys *)
|
|
|
|
(* Calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)
|
|
let rec localisation m_C methode loc_deloc epsilon n prev_critere_D cc=
|
|
|
|
(*Printf.printf "%i\n%!" n;*)
|
|
|
|
(*Util.debug_matrix "m_C" m_C;*)
|
|
|
|
if n == 0
|
|
then m_C
|
|
else
|
|
|
|
(* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)
|
|
let new_m_C m_C methode loc_deloc =
|
|
|
|
(* Fonction de pattern matching en fonction de la méthode *)
|
|
let alphad = m_alpha_d methode m_C
|
|
in
|
|
|
|
(* D critère à maximiser *)
|
|
let critere_D = alphad.d
|
|
in
|
|
|
|
(*Printf.printf "%f\n%!" critere_D;*)
|
|
|
|
(* Matrice des alphas *)
|
|
let m_alpha = alphad.m_alpha
|
|
in
|
|
let norme_alpha = norme m_alpha
|
|
in
|
|
|
|
(*Util.debug_matrix "m_alpha" m_alpha;*)
|
|
|
|
(* alphaij contient le alpha max ainsi que ses indices i et j *)
|
|
let n_rec_alpha = 10 (* Nombre ditération max pour réduire les valeurs de alpha *)
|
|
in
|
|
let alphaij = new_m_alpha m_alpha m_C n_rec_alpha
|
|
in
|
|
|
|
(* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)
|
|
let alpha = (alpha_v loc_deloc alphaij) *. epsilon
|
|
in
|
|
Printf.printf "%i %f %f %f\n%!" n critere_D alpha norme_alpha;
|
|
|
|
|
|
(*Printf.printf "%f\n%!" alpha;*)
|
|
|
|
(* Indice i et j du alpha max après calcul *)
|
|
let indice_i = alphaij.indice_ii
|
|
in
|
|
let indice_j = alphaij.indice_jj
|
|
in
|
|
|
|
(*Printf.printf "%i %i\n%!" indice_i indice_j;*)
|
|
|
|
(* Matrice de rotation *)
|
|
let m_R = f_R alpha
|
|
in
|
|
|
|
(*Util.debug_matrix "m_R" m_R;*)
|
|
|
|
(* Matrice qui va subir la rotation *)
|
|
let m_Ksi = f_Ksi indice_i indice_j m_C
|
|
in
|
|
|
|
(*Util.debug_matrix "m_Ksi" m_Ksi;*)
|
|
|
|
(* Matrice ayant subit la rotation *)
|
|
let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi m_C
|
|
in
|
|
|
|
(*Util.debug_matrix "m_Ksi_tilde" m_Ksi_tilde;*)
|
|
|
|
(* Matrice pour supprimerles coef des orbitales i et j dans la matrice des coef *)
|
|
let m_Psi = f_k m_Ksi indice_i indice_j m_C
|
|
in
|
|
|
|
(*Util.debug_matrix "m_Psi" m_Psi;*)
|
|
|
|
(* Matrice pour ajouter les coef des orbitales i~ et j~ dans la matrice des coef *)
|
|
let m_Psi_tilde = f_k m_Ksi_tilde indice_i indice_j m_C
|
|
in
|
|
|
|
(*Util.debug_matrix "m_Psi_tilde" m_Psi_tilde;*)
|
|
|
|
(* Matrice avec les coef des orbitales i et j remplacés par 0 *)
|
|
let m_interm = f_interm m_C m_Psi
|
|
in
|
|
|
|
(*Util.debug_matrix "m_interm" m_interm;*)
|
|
|
|
(* Matrice après rotation *)
|
|
( Mat.add m_Psi_tilde m_interm, critere_D, norme_alpha, alpha)
|
|
in
|
|
let m_new_m_C , critere_D, norme_alpha, alpha = new_m_C m_C methode loc_deloc
|
|
in
|
|
let _diff = prev_critere_D -. critere_D in
|
|
|
|
(*Util.debug_matrix "new_alpha_m" (f_alpha m_C);*)
|
|
(*Util.debug_matrix "m_new_m_C" m_new_m_C;*)
|
|
|
|
if alpha**2. < cc**2.
|
|
then m_new_m_C
|
|
else
|
|
localisation m_new_m_C methode loc_deloc epsilon (n-1) critere_D cc;;
|
|
|
|
(* Calcul *)
|
|
(* Fonction / Matrice des coef / Méthode("Boys" ou "ER") / Localisation ou non ("loc" ou "deloc") / Pas(<=1.)
|
|
/ Nombre d'itérations max / 0. (valeur de D pour initier la boucle) / critère de convergence sur D*)
|
|
(*
|
|
let new_m_boys = localisation m_C "boys" "loc" 1. 100 0. 10e-7;;
|
|
let new_m_er = localisation m_C "ER" "loc" 1. 100 0. 10e-7;;
|
|
*)
|
|
|
|
|
|
(*Fonction de création d'une list d'entier à partir d'un vecteur de float*)
|
|
let int_list vec =
|
|
let float_list = Vec.to_list vec
|
|
in
|
|
let g a = int_of_float(a)
|
|
in List.map g float_list;;
|
|
|
|
(* Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)
|
|
let miss_elem mat list =
|
|
let n_mo = Mat.dim2 mat
|
|
in
|
|
let vec = Vec.init (n_mo) (fun i ->
|
|
if List.mem i list
|
|
then 0.
|
|
else float_of_int(i))
|
|
in
|
|
let list_int = int_list vec
|
|
in
|
|
List.filter (fun x -> x > 0) list_int;;
|
|
|
|
(* Fonction de séparation d'une matrice en 2 sous matrice, la première matrice correspondant aux colonnes de la matrice dont le numéro est présent
|
|
dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)
|
|
let split_mat mat list =
|
|
let vec_of_mat = Mat.to_col_vecs mat
|
|
in
|
|
let f a = vec_of_mat.(a-1)
|
|
in
|
|
let vec_list_1 = List.map f list
|
|
in
|
|
let list_2 = miss_elem mat list
|
|
in
|
|
let vec_list_2 = List.map f list_2
|
|
in (Mat.of_col_vecs_list vec_list_1,Mat.of_col_vecs_list vec_list_2);;
|
|
|
|
(* Liste des OMs occupées *)
|
|
let list_occ =
|
|
let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))
|
|
in int_list vec_occ;;
|
|
|
|
(* Fonction de rassemblement de 2 matrices *)
|
|
let assemble m_occ m_vir =
|
|
let occ = Mat.to_col_vecs m_occ in
|
|
let vir = Mat.to_col_vecs m_vir in
|
|
Array.concat [ occ ; vir ]
|
|
|> Mat.of_col_vecs
|
|
|
|
(**********************************)
|
|
(*
|
|
m_C;;
|
|
|
|
let list_om = List.init nocc (fun i -> i+1)
|
|
|
|
let m_occ , m_vir = split_mat m_C list_om;;
|
|
|
|
assemble m_occ m_vir;;
|
|
*)
|
|
|
|
|
|
let m_occ , m_vir = split_mat m_C list_occ;;
|
|
|
|
(*
|
|
let new_m_boys = localisation m_C "boys" "loc" 1. 1000 0. 10e-7;;
|
|
let new_m_er = localisation m_C "ER" "loc" 1. 1000 0. 10e-7;;
|
|
*)
|
|
|
|
Printf.printf "Boys\n"
|
|
(*
|
|
let loc_m_occ_boys = localisation m_occ "boys" "loc" 1. 5000 0. 1e-3;;
|
|
let loc_m_vir_boys = localisation m_vir "boys" "loc" 1. 5000 0. 1e-3;;
|
|
let m_assemble_boys = assemble loc_m_occ_boys loc_m_vir_boys;;
|
|
*)
|
|
|
|
(*
|
|
let loc_m_occ_boys = localisation m_occ "boys" "loc" 1. 5000 0. 1e-3;;
|
|
let m_assemble_boys = assemble loc_m_occ_boys m_vir;;
|
|
*)
|
|
|
|
let loc_m_vir_boys = localisation m_vir "boys" "loc" 1. 50000 0. 1e-3;;
|
|
let m_assemble_boys = assemble m_occ loc_m_vir_boys;;
|
|
|
|
Printf.printf "[\n";;
|
|
Mat.as_vec m_assemble_boys
|
|
|> Vec.iter (fun x -> Printf.printf "%20.15e,\n" x);;
|
|
Printf.printf "]\n";;
|
|
|
|
(*
|
|
Printf.printf "ER\n"
|
|
let loc_m_occ_er = localisation m_occ "ER" "loc" 1. 100 0. 1e-3;;
|
|
let loc_m_vir_er = localisation m_vir "ER" "loc" 1. 100 0. 1e-3;;
|
|
let m_assemble_er = assemble loc_m_occ_er loc_m_vir_er;;
|
|
|
|
Mat.as_vec m_assemble_er
|
|
|> Vec.iter (fun x -> Printf.printf "%20.15e\n" x);;
|
|
*)
|
|
|
|
(*
|
|
let mo_base1 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized "Boys")
|
|
~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:new_m_boys ();;
|
|
|
|
let mo_base2 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized "ER")
|
|
~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:new_m_er ();;
|
|
|
|
let mo_base3 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized "Boys")
|
|
~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:m_assemble_boys ();;
|
|
|
|
let mo_base4 = MOBasis.make ~simulation ~mo_type:(MOBasis.Localized "ER")
|
|
~mo_occupation:(MOBasis.mo_occupation mo_basis) ~mo_coef:m_assemble_er ();;
|
|
*)
|
|
|
|
(*
|
|
|
|
(*let mo_basis = MOBasis.of_hartree_fock hf*)
|
|
let mo_basis = mo_base4
|
|
|
|
let ci =
|
|
DeterminantSpace.fci_of_mo_basis mo_basis ~frozen_core:false
|
|
|> CI.make
|
|
|
|
let ci_coef, ci_energy = Lazy.force ci.eigensystem
|
|
let m_ci_coef = Mat.as_vec ci_coef;;
|
|
Vec.iteri (fun i x -> Printf.printf "%d %e\n" i x) m_ci_coef;;
|
|
*)
|
|
|
|
|
|
|
|
|
|
|
|
|