1217 lines
46 KiB
Plaintext
1217 lines
46 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Sauvegarde des anciennes fonctions"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(*\n",
|
|
"\n",
|
|
"(* Localisation de Edminstion ou de Boys *)\n",
|
|
"\n",
|
|
"(* Calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)\n",
|
|
"let rec localisation m_C methode epsilon n prev_critere_D cc=\n",
|
|
"\n",
|
|
" Printf.printf \"%i\\n%!\" n;\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_C\" m_C;*)\n",
|
|
"\n",
|
|
"if n == 0 \n",
|
|
" then m_C\n",
|
|
" else\n",
|
|
" \n",
|
|
" (* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)\n",
|
|
" let new_m_C m_C methode =\n",
|
|
" \n",
|
|
" (* Fonction de pattern matching en fonction de la méthode *)\n",
|
|
" let alphad = m_alpha_d methode m_C \n",
|
|
" in\n",
|
|
" \n",
|
|
" (* D critère à maximiser *)\n",
|
|
" let critere_D = alphad.d \n",
|
|
" in\n",
|
|
" \n",
|
|
" Printf.printf \"%f\\n%!\" critere_D;\n",
|
|
" \n",
|
|
" (* Matrice des alphas *)\n",
|
|
" let m_alpha = alphad.m_alpha\n",
|
|
" in\n",
|
|
" let norme_alpha = norme m_alpha\n",
|
|
" in\n",
|
|
" \n",
|
|
" Printf.printf \"%f\\n%!\" norme_alpha;\n",
|
|
" \n",
|
|
" (*Util.debug_matrix \"m_alpha\" m_alpha;*)\n",
|
|
"\n",
|
|
" (* alphaij contient le alpha max ainsi que ses indices i et j *)\n",
|
|
" let n_rec_alpha = 10 (* Nombre ditération max pour réduire les valeurs de alpha *)\n",
|
|
" in\n",
|
|
" let alphaij = new_m_alpha m_alpha m_C n_rec_alpha \n",
|
|
" in\n",
|
|
"\n",
|
|
" (* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)\n",
|
|
" let alpha = (alphaij.alpha_max) *. epsilon \n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Printf.printf \"%f\\n%!\" alpha;*)\n",
|
|
" \n",
|
|
" (* Indice i et j du alpha max après calcul *)\n",
|
|
" let indice_i = alphaij.indice_ii \n",
|
|
" in\n",
|
|
" let indice_j = alphaij.indice_jj \n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Printf.printf \"%i %i\\n%!\" indice_i indice_j;*)\n",
|
|
" \n",
|
|
" (* Matrice de rotation *)\n",
|
|
" let m_R = f_R alpha \n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_R\" m_R;*)\n",
|
|
"\n",
|
|
" (* Matrice qui va subir la rotation *)\n",
|
|
" let m_Ksi = f_Ksi indice_i indice_j m_C \n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_Ksi\" m_Ksi;*)\n",
|
|
"\n",
|
|
" (* Matrice ayant subit la rotation *)\n",
|
|
" let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi m_C \n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;*)\n",
|
|
"\n",
|
|
" (* Matrice pour supprimerles coef des orbitales i et j dans la matrice des coef *)\n",
|
|
" let m_Psi = f_k m_Ksi indice_i indice_j m_C\n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_Psi\" m_Psi;*)\n",
|
|
"\n",
|
|
" (* Matrice pour ajouter les coef des orbitales i~ et j~ dans la matrice des coef *)\n",
|
|
" let m_Psi_tilde = f_k m_Ksi_tilde indice_i indice_j m_C \n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_Psi_tilde\" m_Psi_tilde;*)\n",
|
|
"\n",
|
|
" (* Matrice avec les coef des orbitales i et j remplacés par 0 *)\n",
|
|
" let m_interm = f_interm m_C m_Psi \n",
|
|
" in\n",
|
|
" \n",
|
|
" (*Util.debug_matrix \"m_interm\" m_interm;*)\n",
|
|
" \n",
|
|
" (* Matrice après rotation *)\n",
|
|
" ( Mat.add m_Psi_tilde m_interm, critere_D)\n",
|
|
" \n",
|
|
" in\n",
|
|
" let m_new_m_C , critere_D = new_m_C m_C methode \n",
|
|
" in\n",
|
|
" let diff = prev_critere_D -. critere_D\n",
|
|
"\n",
|
|
"in\n",
|
|
"\n",
|
|
"(*Util.debug_matrix \"new_alpha_m\" (f_alpha m_C);*)\n",
|
|
"(*Util.debug_matrix \"m_new_m_C\" m_new_m_C;*)\n",
|
|
"\n",
|
|
"if diff**2. < cc**2.\n",
|
|
" then m_new_m_C\n",
|
|
" else\n",
|
|
"localisation m_new_m_C methode epsilon (n-1) critere_D cc;;\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(* Fonction de calcul de D Boys *)\n",
|
|
"(*\n",
|
|
"let d_boys m_C = \n",
|
|
"\n",
|
|
" let phi_x_phi =\n",
|
|
" Multipole.matrix_x multipoles \n",
|
|
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C \n",
|
|
" in\n",
|
|
" \n",
|
|
" (*Util.debug_matrix \"phi_x_phi\" phi_x_phi;*)\n",
|
|
" \n",
|
|
" let phi_y_phi =\n",
|
|
" Multipole.matrix_y multipoles \n",
|
|
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C \n",
|
|
" in\n",
|
|
" \n",
|
|
" (*Util.debug_matrix \"phi_y_phi\" phi_y_phi;*)\n",
|
|
" \n",
|
|
" let phi_z_phi =\n",
|
|
" Multipole.matrix_z multipoles \n",
|
|
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
|
|
"\n",
|
|
" in\n",
|
|
" \n",
|
|
" (*Util.debug_matrix \"phi_z_phi\" phi_z_phi;*)\n",
|
|
" \n",
|
|
" let v_D_boys = \n",
|
|
" let n_mo = Mat.dim2 m_C\n",
|
|
" in\n",
|
|
" Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)\n",
|
|
"in\n",
|
|
"Vec.sum v_D_boys;;\n",
|
|
"*)\n",
|
|
"\n",
|
|
"(*************************)\n",
|
|
"(*\n",
|
|
"Multipole.matrix_x multipoles;;\n",
|
|
"d_boys m_C;;\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(* Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)\n",
|
|
"(*\n",
|
|
"let miss_elem mat list = \n",
|
|
" let n_mo = Mat.dim2 mat\n",
|
|
" in\n",
|
|
" let vec = Vec.init (n_mo) (fun i ->\n",
|
|
" if List.mem i list\n",
|
|
" then 0.\n",
|
|
" else float_of_int(i))\n",
|
|
" in\n",
|
|
" let vec_list = Vec.to_list vec\n",
|
|
" in\n",
|
|
" let g a = int_of_float(a)\n",
|
|
" in\n",
|
|
" let vec_list_int = List.map g vec_list\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n",
|
|
"(*\n",
|
|
"let f_alpha m_C =\n",
|
|
"\n",
|
|
" let n_mo = Mat.dim2 m_C in\n",
|
|
" let t0 = Sys.time () in\n",
|
|
" \n",
|
|
" let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
|
|
" let m_a12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
|
|
" \n",
|
|
" (* Tableaux temporaires *)\n",
|
|
" let m_pqr =\n",
|
|
" Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)\n",
|
|
" in\n",
|
|
" let m_qr_i = Mat.create (n_ao*n_ao) n_mo in\n",
|
|
" let m_ri_j = Mat.create (n_ao*n_mo) n_mo in\n",
|
|
" let m_ij_k = Mat.create (n_mo*n_mo) n_mo in\n",
|
|
" \n",
|
|
" Array.iter (fun s ->\n",
|
|
" (* Grosse boucle externe sur s *)\n",
|
|
" Array.iter (fun r ->\n",
|
|
" Array.iter (fun q ->\n",
|
|
" Array.iter (fun p ->\n",
|
|
" m_pqr.{p,q,r} <- ERI.get_phys ee_ints p q r s\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" ) (Util.array_range 1 n_ao);\n",
|
|
" \n",
|
|
" (* Conversion d'un tableau a 3 indices en une matrice nao x nao^2 *)\n",
|
|
" let m_p_qr =\n",
|
|
" Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]\n",
|
|
" |> Bigarray.array2_of_genarray\n",
|
|
" in\n",
|
|
" \n",
|
|
" let m_qr_i =\n",
|
|
" (* (qr,i) = <i r|q s> = \\sum_p <p r | q s> C_{pi} *)\n",
|
|
" gemm ~transa:`T ~c:m_qr_i m_p_qr m_C\n",
|
|
" in\n",
|
|
" \n",
|
|
" let m_q_ri =\n",
|
|
" (* Transformation de la matrice (qr,i) en (q,ri) *)\n",
|
|
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo)\n",
|
|
" in\n",
|
|
" \n",
|
|
" let m_ri_j =\n",
|
|
" (* (ri,j) = <i r | j s> = \\sum_q <i r | q s> C_{bj} *)\n",
|
|
" gemm ~transa:`T ~c:m_ri_j m_q_ri m_C\n",
|
|
" in\n",
|
|
" \n",
|
|
" let m_r_ij =\n",
|
|
" (* Transformation de la matrice (ri,j) en (r,ij) *)\n",
|
|
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo)\n",
|
|
" in\n",
|
|
" \n",
|
|
" let m_ij_k =\n",
|
|
" (* (ij,k) = <i k | j s> = \\sum_r <i r | j s> C_{rk} *)\n",
|
|
" gemm ~transa:`T ~c:m_ij_k m_r_ij m_C\n",
|
|
" in\n",
|
|
" \n",
|
|
" let m_ijk =\n",
|
|
" (* Transformation de la matrice (ei,j) en (e,ij) *)\n",
|
|
" Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |]\n",
|
|
" |> Bigarray.array3_of_genarray\n",
|
|
" in\n",
|
|
" \n",
|
|
" Array.iter (fun j ->\n",
|
|
" Array.iter (fun i ->\n",
|
|
" m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});\n",
|
|
" m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j} -.\n",
|
|
" 0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_C.{s,i} +.\n",
|
|
" (m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_C.{s,j})\n",
|
|
" ) (Util.array_range 1 n_mo)\n",
|
|
" ) (Util.array_range 1 n_mo)\n",
|
|
" ) (Util.array_range 1 n_ao);\n",
|
|
" \n",
|
|
" let t1 = Sys.time () in\n",
|
|
" Printf.printf \"t = %f s\\n%!\" (t1 -. t0);\n",
|
|
" Mat.init_cols n_mo n_mo ( fun i j ->\n",
|
|
" if i= j then 0.\n",
|
|
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2. ))))\n",
|
|
" );;\n",
|
|
"\n",
|
|
"(*********************)\n",
|
|
"\n",
|
|
"f_alpha m_C;;\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(*\n",
|
|
"(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n",
|
|
"let f_alpha m_C eps=\n",
|
|
" let n_mo = Mat.dim2 m_C\n",
|
|
" in\n",
|
|
" (* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)\n",
|
|
" let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
|
|
" let m_a12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
|
|
" Array.iter (fun a ->\n",
|
|
" Array.iter (fun b ->\n",
|
|
" let mca = Vec.init n_mo (fun i -> m_C.{a,i} *. m_C.{b,i}) in\n",
|
|
" Array.iter (fun e ->\n",
|
|
" Array.iter (fun f ->\n",
|
|
" let integral = ERI.get_phys ee_ints f b e a in\n",
|
|
" if abs_float integral > eps then\n",
|
|
" Array.iter ( fun i -> \n",
|
|
" let mcei = m_C.{e,i} *. m_C.{f,i} in\n",
|
|
" Array.iter ( fun j -> \n",
|
|
" let x = m_C.{e,i} *. m_C.{f,j} *. integral in\n",
|
|
" let mcaij = ( mca.{i} -. mca.{j} ) in\n",
|
|
" m_b12.{i,j} <- m_b12.{i,j} +. mcaij *. x;\n",
|
|
" m_a12.{i,j} <- m_a12.{i,j} +. m_C.{a,i} *. m_C.{b,j} *. x\n",
|
|
" -. 0.25 *. ( mcei -. m_C.{e,j} *. m_C.{f,j} ) *. mcaij *. integral \n",
|
|
" ) (Util.array_range 1 n_mo)\n",
|
|
" ) (Util.array_range 1 n_mo)\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" ) (Util.array_range 1 n_ao);\n",
|
|
"\n",
|
|
" Mat.init_cols n_mo n_mo ( fun i j -> \n",
|
|
" if i= j then 0. \n",
|
|
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))));;\n",
|
|
"\n",
|
|
"(*********************)\n",
|
|
"\n",
|
|
"f_alpha m_C 1.e-5;; \n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(*\n",
|
|
"(*Fonction général de calcul des intégrales*) \n",
|
|
"let integral_general g i j =\n",
|
|
"Array.map (fun a ->\n",
|
|
" let v = \n",
|
|
" Array.map (fun b ->\n",
|
|
" let u = \n",
|
|
" Array.map (fun e ->\n",
|
|
" let t = Array.map (fun f ->\n",
|
|
" (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum t\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum u\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum v\n",
|
|
") (Util.array_range 1 n_ao)\n",
|
|
"|> sum \n",
|
|
"\n",
|
|
"\n",
|
|
"(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n",
|
|
"\n",
|
|
"let f_alpha m_C =\n",
|
|
" let n_mo = Mat.dim2 m_C\n",
|
|
" in\n",
|
|
" (* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)\n",
|
|
" let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> \n",
|
|
" integral_general (fun a b e f i j ->\n",
|
|
" ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}\n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
"\n",
|
|
" in\n",
|
|
" (* Fonction de calcul de toutes les intégrales A_12 -> Matrice, dépend de m_C *)\n",
|
|
" let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->\n",
|
|
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j} *. m_C.{e,i} *. m_C.{f,j} \n",
|
|
" -. 0.25 *. (( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} ) \n",
|
|
" *. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ))\n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
"in\n",
|
|
"Mat.init_cols n_mo n_mo ( fun i j -> \n",
|
|
" if i= j \n",
|
|
" then 0. \n",
|
|
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))));;\n",
|
|
"\n",
|
|
"(*********************)\n",
|
|
"(*\n",
|
|
"f_alpha m_C;; \n",
|
|
"*)\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(*\n",
|
|
"(* Calcul de D -> critère à maximiser dans ER*)\n",
|
|
"let s_D m_C = \n",
|
|
" let v_D = \n",
|
|
" let n_mo = Mat.dim2 m_C\n",
|
|
" in\n",
|
|
" let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
|
|
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
" in Vec.init n_mo ( fun i -> m_D.{i,i} )\n",
|
|
" in Vec.sum v_D ;;\n",
|
|
"\n",
|
|
"(******************)\n",
|
|
"let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
|
|
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
|
|
" ) i j\n",
|
|
" );;\n",
|
|
"let toto = s_D m_C;;\n",
|
|
"\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Boucle localisation avec fonctions définies à l'extérieur"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(* Ancienne boucle\n",
|
|
"\n",
|
|
"(* Localisation de Edminstion *)\n",
|
|
"let n_rec_alpha = 10;;\n",
|
|
"(* Calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)\n",
|
|
"let rec final_m_C m_C n =\n",
|
|
"\n",
|
|
" Printf.printf \"%i\\n%!\" n;\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_C\" m_C;*)\n",
|
|
"\n",
|
|
"if n == 0 \n",
|
|
" then m_C\n",
|
|
" else\n",
|
|
" \n",
|
|
" (* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)\n",
|
|
" let new_m_C m_C =\n",
|
|
" \n",
|
|
" (* D critère à maximiser *)\n",
|
|
" let critere_D = s_D m_C (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" Printf.printf \"%f\\n%!\" critere_D;\n",
|
|
" \n",
|
|
" (* Matrice des alphas *)\n",
|
|
" let m_alpha = f_alpha m_C (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_alpha\" m_alpha;*)\n",
|
|
"\n",
|
|
" (* alphaij contient le alpha max ainsi que ses indices i et j *)\n",
|
|
" let alphaij = new_m_alpha m_alpha n_rec_alpha (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
"\n",
|
|
" (* Valeur de alpha max après calcul *)\n",
|
|
" let alpha = alphaij.alpha_max (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
"\n",
|
|
" (* Indice i et j du alpha max après calcul *)\n",
|
|
" let indice_i = alphaij.indice_ii (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" let indice_j = alphaij.indice_jj (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Printf.printf \"%i %i\\n%!\" indice_i indice_j;*)\n",
|
|
"\n",
|
|
" (* Matrice de rotaion *)\n",
|
|
" let m_R = f_R alpha (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_R\" m_R;*)\n",
|
|
"\n",
|
|
" (* Matrice qui va subir la rotation *)\n",
|
|
" let m_Ksi = f_Ksi indice_i indice_j m_C (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_Ksi\" m_Ksi;*)\n",
|
|
"\n",
|
|
" (* Matrice ayant subit la rotation *)\n",
|
|
" let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;*)\n",
|
|
"\n",
|
|
" (* Matrice pour supprimerles coef des orbitales i et j dans la matrice des coef *)\n",
|
|
" let m_Psi = f_Psi m_Ksi indice_i indice_j (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_Psi\" m_Psi;*)\n",
|
|
"\n",
|
|
" (* Matrice pour ajouter les coef des orbitales i~ et j~ dans la matrice des coef *)\n",
|
|
" let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
"\n",
|
|
" (*Util.debug_matrix \"m_Psi_tilde\" m_Psi_tilde;*)\n",
|
|
"\n",
|
|
" (* Matrice avec les coef des orbitales i et j remplacés par 0 *)\n",
|
|
" let m_interm = f_interm m_C m_Psi (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" (* Matrice après rotation *)\n",
|
|
" Mat.add m_Psi_tilde m_interm\n",
|
|
" in\n",
|
|
" let m_new_m_C = new_m_C m_C (* Fonction -> constante *)\n",
|
|
"\n",
|
|
"in\n",
|
|
"\n",
|
|
"(*Util.debug_matrix \"new_alpha_m\" (f_alpha m_C);*)\n",
|
|
"(*Util.debug_matrix \"m_new_m_C\" m_new_m_C;*)\n",
|
|
"\n",
|
|
"final_m_C m_new_m_C (n-1);;\n",
|
|
"\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(*\n",
|
|
"(*Cellule de calcul*)\n",
|
|
"final_m_C m_C 10;;\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(*Définition des fonctions teste*)\n",
|
|
"(*\n",
|
|
"(* Localisation de Edminstion *)\n",
|
|
"\n",
|
|
"(* Définitions de base nécessaire pour la suite *)\n",
|
|
"let ee_ints = AOBasis.ee_ints ao_basis;;\n",
|
|
"let m_C = MOBasis.mo_coef mo_basis;;\n",
|
|
"let n_ao = Mat.dim1 m_C ;;\n",
|
|
"let n_mo = Mat.dim2 m_C ;;\n",
|
|
"let sum a = \n",
|
|
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
|
|
" \n",
|
|
"(******************************************************************************************************)\n",
|
|
"\n",
|
|
"(*Fonction général de calcul des intégrales*) \n",
|
|
"let integral_general g i j =\n",
|
|
"Array.map (fun a ->\n",
|
|
" let v = \n",
|
|
" Array.map (fun b ->\n",
|
|
" let u = \n",
|
|
" Array.map (fun e ->\n",
|
|
" let t = Array.map (fun f ->\n",
|
|
" (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum t\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum u\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum v\n",
|
|
") (Util.array_range 1 n_ao)\n",
|
|
"|> sum \n",
|
|
";;\n",
|
|
"\n",
|
|
"(******************************************************************************************************)\n",
|
|
"\n",
|
|
"(* Fonction de calcul de tous les alpha -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n",
|
|
"let f_alpha m_C =\n",
|
|
"\n",
|
|
"(* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)\n",
|
|
"let m_b12 = Mat.init_cols n_ao n_ao (fun i j -> \n",
|
|
"integral_general (fun a b e f i j ->\n",
|
|
" ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}\n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
"\n",
|
|
"in\n",
|
|
"(* Fonction de calcul de toutes les intégrales A_12 -> Matrice, dépend de m_C *)\n",
|
|
"let m_a12 = Mat.init_cols n_ao n_ao (fun i j ->\n",
|
|
"integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j} *. m_C.{e,i} *. m_C.{f,j} \n",
|
|
" -. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} ) \n",
|
|
" *. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )\n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
"in\n",
|
|
" Mat.init_cols n_ao n_ao ( fun i j ->\n",
|
|
" if i= j \n",
|
|
" then 0. \n",
|
|
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))))\n",
|
|
"\n",
|
|
";;\n",
|
|
"\n",
|
|
"\n",
|
|
"f_alpha m_C;;\n",
|
|
" \n",
|
|
"(******************************************************************************************************)\n",
|
|
"\n",
|
|
"(* Calcul de D *)\n",
|
|
" let s_D m_C = \n",
|
|
" let v_D = \n",
|
|
" let m_D = Mat.init_cols n_ao n_ao (fun i j ->\n",
|
|
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
" in Vec.init n_ao ( fun i -> m_D.{i,1} )\n",
|
|
" in Vec.sum v_D ;;\n",
|
|
" \n",
|
|
" s_D m_C;;\n",
|
|
" \n",
|
|
"(******************************************************************************************************)\n",
|
|
"\n",
|
|
"(* Si alpha max > pi/2 on doit soustraire pi/2 à la matrice des alphas *)\n",
|
|
"let m_alpha = f_alpha m_C;;\n",
|
|
"\n",
|
|
"type alphaij = {\n",
|
|
"alpha_max : float;\n",
|
|
"indice_ii : int;\n",
|
|
"indice_jj : int;}\n",
|
|
"\n",
|
|
"let rec new_m_alpha m_alpha n_rec_alpha=\n",
|
|
" let alpha_m =\n",
|
|
" Printf.printf \"%i\\n%!\" n_rec_alpha;\n",
|
|
" if n_rec_alpha == 0 \n",
|
|
" then m_alpha \n",
|
|
" else Mat.init_cols n_ao n_ao (fun i j -> \n",
|
|
" if (m_alpha.{i,j}) > 3.14 /. 2. \n",
|
|
" then (m_alpha.{i,j} -. ( 3.14 /. 2.))\n",
|
|
" else if m_alpha.{i,j} < -. 3.14 /. 2.\n",
|
|
" then (m_alpha.{i,j} +. ( 3.14 /. 2.))\n",
|
|
" else if m_alpha.{i,j} < 0. \n",
|
|
" then -. m_alpha.{i,j}\n",
|
|
" else m_alpha.{i,j} )\n",
|
|
" \n",
|
|
" in \n",
|
|
" Util.debug_matrix \"alpha_m\" alpha_m;\n",
|
|
" let max_element3 alpha_m = \n",
|
|
" Mat.as_vec alpha_m\n",
|
|
" |> iamax\n",
|
|
" in\n",
|
|
" let indice_ii = \n",
|
|
" let max = max_element3 alpha_m\n",
|
|
" in\n",
|
|
" Printf.printf \"%i\\n%!\" max;\n",
|
|
" (max - 1) mod n_ao +1 \n",
|
|
" in\n",
|
|
" let indice_jj = \n",
|
|
" let max = max_element3 alpha_m\n",
|
|
" in\n",
|
|
" (max - 1) / n_ao +1\n",
|
|
" in\n",
|
|
" let alpha alpha_m = \n",
|
|
" let i = indice_ii \n",
|
|
" in\n",
|
|
" let j = indice_jj \n",
|
|
" in\n",
|
|
" alpha_m.{i,j}\n",
|
|
" in\n",
|
|
" let alpha_max = alpha alpha_m \n",
|
|
" in\n",
|
|
" Printf.printf \"%f\\n%!\" alpha_max;\n",
|
|
" if alpha_max < 3.14 /. 2.\n",
|
|
" then {alpha_max; indice_ii; indice_jj}\n",
|
|
" else new_m_alpha alpha_m (n_rec_alpha-1);;\n",
|
|
" \n",
|
|
"\n",
|
|
"new_m_alpha m_alpha 10 ;;\n",
|
|
"\n",
|
|
"let alphaij = new_m_alpha m_alpha 10 ;;\n",
|
|
"\n",
|
|
"(******************************************************************************************)\n",
|
|
"let alpha = alphaij.alpha_max\n",
|
|
"\n",
|
|
"(* Matrice de rotation 2 par 2 *)\n",
|
|
"let m_R alpha =\n",
|
|
" Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha\n",
|
|
" else if i>j then sin alpha \n",
|
|
" else -. sin alpha);;\n",
|
|
"\n",
|
|
"\n",
|
|
"(******************************************************************************************)\n",
|
|
"\n",
|
|
"let indice_i = alphaij.indice_ii;;\n",
|
|
"let indice_j = alphaij.indice_jj;;\n",
|
|
"let m_R = m_R alpha;;\n",
|
|
"\n",
|
|
"(* 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)\n",
|
|
"pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n",
|
|
"let m_Ksi indice_i indice_j m_C = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_C.{i,indice_i} else m_C.{i,indice_j} );;\n",
|
|
"\n",
|
|
"m_Ksi indice_i indice_j m_C;;\n",
|
|
"\n",
|
|
"let m_Ksi = m_Ksi indice_i indice_j m_C;;\n",
|
|
"\n",
|
|
"(* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle\n",
|
|
"on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n",
|
|
"let m_Ksi_tilde m_R m_Ksi = gemm m_Ksi m_R;;\n",
|
|
"\n",
|
|
"m_Ksi_tilde m_R m_Ksi ;;\n",
|
|
"\n",
|
|
"(******************************************************************************************)\n",
|
|
"let m_Ksi_tilde = m_Ksi_tilde m_R m_Ksi\n",
|
|
"\n",
|
|
"(* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur \n",
|
|
"les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice\n",
|
|
"des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec \n",
|
|
"celle contenant i~ et j~ *)\n",
|
|
"\n",
|
|
"(* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n",
|
|
"let m_Psi_tilde m_Ksi_tilde indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi_tilde.{i,1}\n",
|
|
" else if j=indice_j then m_Ksi_tilde.{i,2}\n",
|
|
" else 0.);;\n",
|
|
"\n",
|
|
"m_Psi_tilde m_Ksi_tilde indice_i indice_j;;\n",
|
|
"\n",
|
|
"(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
|
|
"let m_Psi m_Ksi indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi.{i,1}\n",
|
|
" else if j=indice_j then m_Ksi.{i,2}\n",
|
|
" else 0.);;\n",
|
|
" \n",
|
|
"m_Psi m_Ksi indice_i indice_j;;\n",
|
|
"\n",
|
|
"(******************************************************************************************)\n",
|
|
"let m_Psi = m_Psi m_Ksi indice_i indice_j;;\n",
|
|
"let m_Psi_tilde = m_Psi_tilde m_Ksi_tilde indice_i indice_j;; \n",
|
|
"\n",
|
|
"(* 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\n",
|
|
"par la matrice *)\n",
|
|
"let m_interm m_C m_Psi = Mat.sub m_C m_Psi;;\n",
|
|
"\n",
|
|
"let m_interm = m_interm m_C m_Psi;;\n",
|
|
"\n",
|
|
"(* Construction de la nouvelle matrice d'OMs phi~ *)\n",
|
|
"let m_Phi_tilde m_Psi_tilde m_interm = Mat.add m_Psi_tilde m_interm;;\n",
|
|
"\n",
|
|
"m_Phi_tilde m_Psi_tilde m_interm;;\n",
|
|
"\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Programmes avec fonctions définies à l'intérieurs"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(* Localisation de Edminstion *)\n",
|
|
"\n",
|
|
"(*Fonctionne -> programme avec les fonctions définies à l'intérieur*)\n",
|
|
"\n",
|
|
"(*\n",
|
|
"\n",
|
|
"(* Définitions de base nécessaire pour la suite *)\n",
|
|
"let ee_ints = AOBasis.ee_ints ao_basis;;\n",
|
|
"let m_C = MOBasis.mo_coef mo_basis;;\n",
|
|
"let n_ao = Mat.dim1 m_C ;;\n",
|
|
"let n_mo = Mat.dim2 m_C ;;\n",
|
|
"let sum a = \n",
|
|
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
|
|
" \n",
|
|
"(******************************************************************************************************)\n",
|
|
"\n",
|
|
"(*Fonction général de calcul des intégrales*) \n",
|
|
"let integral_general g i j =\n",
|
|
"Array.map (fun a ->\n",
|
|
" let v = \n",
|
|
" Array.map (fun b ->\n",
|
|
" let u = \n",
|
|
" Array.map (fun e ->\n",
|
|
" let t = Array.map (fun f ->\n",
|
|
" (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum t\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum u\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum v\n",
|
|
") (Util.array_range 1 n_ao)\n",
|
|
"|> sum \n",
|
|
";;\n",
|
|
"\n",
|
|
"\n",
|
|
"type alphaij = {\n",
|
|
" alpha_max : float;\n",
|
|
" indice_ii : int;\n",
|
|
" indice_jj : int;};;\n",
|
|
" \n",
|
|
"let n_rec_alpha = 10;;\n",
|
|
"\n",
|
|
"(*\n",
|
|
"let alpha = 1.;;\n",
|
|
"let y_R alpha =\n",
|
|
" Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha\n",
|
|
" else if i>j then sin alpha \n",
|
|
" else -. sin alpha);;\n",
|
|
"let y_R = y_R alpha;; \n",
|
|
"let m_C = gemm m_C y_R;;\n",
|
|
"*)\n",
|
|
"(******************************************************************************************************)\n",
|
|
"\n",
|
|
"let rec final_m_C m_C n =\n",
|
|
"\n",
|
|
"Printf.printf \"%i\\n%!\" n;\n",
|
|
"(*\n",
|
|
"Util.debug_matrix \"m_C\" m_C;\n",
|
|
"*)\n",
|
|
"\n",
|
|
"if n == 0 then m_C\n",
|
|
" else\n",
|
|
"\n",
|
|
"let new_m_C m_C =\n",
|
|
"\n",
|
|
" (* Fonction de calcul de tous les alpha -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n",
|
|
" let f_alpha m_C =\n",
|
|
"\n",
|
|
" (* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)\n",
|
|
" let m_b12 = Mat.init_cols n_ao n_ao (fun i j -> \n",
|
|
" integral_general (fun a b e f i j ->\n",
|
|
" ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}\n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
"\n",
|
|
" in\n",
|
|
" (* Fonction de calcul de toutes les intégrales A_12 -> Matrice, dépend de m_C *)\n",
|
|
" let m_a12 = Mat.init_cols n_ao n_ao (fun i j ->\n",
|
|
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j} *. m_C.{e,i} *. m_C.{f,j} \n",
|
|
" -. 0.25 *. (( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} ) \n",
|
|
" *. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ))\n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
" in\n",
|
|
" Mat.init_cols n_ao n_ao ( fun i j -> if i= j \n",
|
|
" then 0. \n",
|
|
" else\n",
|
|
" 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))))\n",
|
|
" in\n",
|
|
" \n",
|
|
" (* Calcul de D *)\n",
|
|
" let s_D m_C = \n",
|
|
" let v_D = \n",
|
|
" let m_D = Mat.init_cols n_ao n_ao (fun i j ->\n",
|
|
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
" in Vec.init n_ao ( fun i -> m_D.{i,1} )\n",
|
|
" in Vec.sum v_D \n",
|
|
" in\n",
|
|
" let critere_D = s_D m_C (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" Printf.printf \"%f\\n%!\" critere_D;\n",
|
|
" \n",
|
|
" let m_alpha = f_alpha m_C (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" (*\n",
|
|
" Util.debug_matrix \"m_alpha\" m_alpha;\n",
|
|
" *)\n",
|
|
" \n",
|
|
" (* Détermination alpha_max et ses indices i et j.\n",
|
|
" Si alpha max > pi/2 on soustrait pi/2 à la matrice des alphas de manière récursive *)\n",
|
|
" let rec new_m_alpha m_alpha n_rec_alpha=\n",
|
|
" let alpha_m =\n",
|
|
" \n",
|
|
" Printf.printf \"%i\\n%!\" n_rec_alpha;\n",
|
|
" \n",
|
|
" if n_rec_alpha == 0 \n",
|
|
" then m_alpha \n",
|
|
" else Mat.init_cols n_ao n_ao (fun i j -> \n",
|
|
" if (m_alpha.{i,j}) > (3.14 /. 2.) \n",
|
|
" then (m_alpha.{i,j} -. ( 3.14 /. 2.))\n",
|
|
" else if m_alpha.{i,j} < -. 3.14 /. 2.\n",
|
|
" then (m_alpha.{i,j} +. ( 3.14 /. 2.))\n",
|
|
" else if m_alpha.{i,j} < 0. \n",
|
|
" then -. m_alpha.{i,j}\n",
|
|
" else m_alpha.{i,j} )\n",
|
|
" in \n",
|
|
" \n",
|
|
" Util.debug_matrix \"alpha_m\" alpha_m;\n",
|
|
" \n",
|
|
" (* Détermination de l'emplacement du alpha max *)\n",
|
|
" let max_element3 alpha_m = \n",
|
|
" Mat.as_vec alpha_m\n",
|
|
" |> iamax\n",
|
|
" in\n",
|
|
" \n",
|
|
" (* indice i du alpha max *)\n",
|
|
" let indice_ii = \n",
|
|
" let max = max_element3 alpha_m (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" (*\n",
|
|
" Printf.printf \"%i\\n%!\" max;\n",
|
|
" *)\n",
|
|
" (max - 1) mod n_ao +1 \n",
|
|
" in\n",
|
|
" \n",
|
|
" (* indice j du alpha max *)\n",
|
|
" let indice_jj = \n",
|
|
" let max = max_element3 alpha_m (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" (max - 1) / n_ao +1\n",
|
|
" in\n",
|
|
" \n",
|
|
" (* Valeur du alpha max*)\n",
|
|
" let alpha alpha_m = \n",
|
|
" let i = indice_ii \n",
|
|
" in\n",
|
|
" let j = indice_jj \n",
|
|
" in\n",
|
|
" (*\n",
|
|
" Printf.printf \"%i %i\\n%!\" i j;\n",
|
|
" *)\n",
|
|
" alpha_m.{i,j}\n",
|
|
" \n",
|
|
" in\n",
|
|
" let alpha_max = alpha alpha_m (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" Printf.printf \"%f\\n%!\" alpha_max;\n",
|
|
" if alpha_max < 3.14 /. 2.\n",
|
|
" then {alpha_max; indice_ii; indice_jj}\n",
|
|
" else new_m_alpha alpha_m (n_rec_alpha-1)\n",
|
|
" in\n",
|
|
" let alphaij = new_m_alpha m_alpha n_rec_alpha (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" \n",
|
|
" (* Valeur de alpha max après calcul *)\n",
|
|
" let alpha = alphaij.alpha_max (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" \n",
|
|
" (* Matrice de rotation 2 par 2 *)\n",
|
|
" let f_R alpha =\n",
|
|
" Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha\n",
|
|
" else if i>j then sin alpha \n",
|
|
" else -. sin alpha )\n",
|
|
" in\n",
|
|
" (* Indice i et j du alpha max après calcul *)\n",
|
|
" let indice_i = alphaij.indice_ii (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" let indice_j = alphaij.indice_jj (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" \n",
|
|
" Printf.printf \"%i %i\\n%!\" indice_i indice_j;\n",
|
|
" \n",
|
|
" let m_R = f_R alpha (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" Util.debug_matrix \"m_R\" m_R;\n",
|
|
" (* 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)\n",
|
|
" pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n",
|
|
" let f_Ksi indice_i indice_j m_C = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_C.{i,indice_i} else m_C.{i,indice_j} )\n",
|
|
" in\n",
|
|
" let m_Ksi = f_Ksi indice_i indice_j m_C (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" \n",
|
|
" Util.debug_matrix \"m_Ksi\" m_Ksi;\n",
|
|
" \n",
|
|
" (* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle\n",
|
|
" on obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n",
|
|
" let f_Ksi_tilde m_R m_Ksi = gemm m_Ksi m_R\n",
|
|
" in\n",
|
|
" let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" \n",
|
|
" Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n",
|
|
" \n",
|
|
" \n",
|
|
" (* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur \n",
|
|
" les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice\n",
|
|
" des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec \n",
|
|
" celle contenant i~ et j~ *)\n",
|
|
" \n",
|
|
" (* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n",
|
|
" let f_Psi_tilde m_Ksi_tilde indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi_tilde.{i,1}\n",
|
|
" else if j=indice_j then m_Ksi_tilde.{i,2}\n",
|
|
" else 0.)\n",
|
|
"\n",
|
|
" in\n",
|
|
"\n",
|
|
" (* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
|
|
" let f_Psi m_Ksi indice_i indice_j = Mat.init_cols n_ao n_ao (fun i j -> if j=indice_i then m_Ksi.{i,1}\n",
|
|
" else if j=indice_j then m_Ksi.{i,2}\n",
|
|
" else 0.)\n",
|
|
" in\n",
|
|
" let m_Psi = f_Psi m_Ksi indice_i indice_j (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" \n",
|
|
" Util.debug_matrix \"m_Psi\" m_Psi;\n",
|
|
" \n",
|
|
" let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j (* Fonction -> constante *)\n",
|
|
" in\n",
|
|
" \n",
|
|
" Util.debug_matrix \"m_Psi_tilde\" m_Psi_tilde;\n",
|
|
" \n",
|
|
" (* 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\n",
|
|
" par la matrice *)\n",
|
|
" let f_interm m_C m_Psi = Mat.sub m_C m_Psi\n",
|
|
" in\n",
|
|
" let m_interm = f_interm m_C m_Psi (* Fonction -> constante *)\n",
|
|
"in\n",
|
|
"Mat.add m_Psi_tilde m_interm\n",
|
|
"in\n",
|
|
"let m_new_m_C = new_m_C m_C (* Fonction -> constante *)\n",
|
|
"in\n",
|
|
"Util.debug_matrix \"new_alpha_m\" (f_alpha m_C);\n",
|
|
"\n",
|
|
"Util.debug_matrix \"m_new_m_C\" m_new_m_C;\n",
|
|
"\n",
|
|
"final_m_C m_new_m_C (n-1);;\n",
|
|
"\n",
|
|
"(*****************************)\n",
|
|
"\n",
|
|
"final_m_C m_C 10;;\n",
|
|
"\n",
|
|
"*)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"(*\n",
|
|
"\n",
|
|
"(* Localisation de Edminstion *)\n",
|
|
"let ee_ints = AOBasis.ee_ints ao_basis;;\n",
|
|
"let m_C = MOBasis.mo_coef mo_basis;;\n",
|
|
"let n_ao = Mat.dim1 m_C ;;\n",
|
|
"let n_mo = Mat.dim2 m_C ;;\n",
|
|
"let sum a = \n",
|
|
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
|
|
";;\n",
|
|
"\n",
|
|
"(*Calcul des intégrales*) \n",
|
|
"let integral_general g i j =\n",
|
|
"Array.map (fun a ->\n",
|
|
" let v = \n",
|
|
" Array.map (fun b ->\n",
|
|
" let u = \n",
|
|
" Array.map (fun e ->\n",
|
|
" let t = Array.map (fun f ->\n",
|
|
" (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum t\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum u\n",
|
|
" ) (Util.array_range 1 n_ao)\n",
|
|
" in sum v\n",
|
|
") (Util.array_range 1 n_ao)\n",
|
|
"|> sum\n",
|
|
";;\n",
|
|
"\n",
|
|
"(* Calcul de toutes les intégrales B_12 -> Matrice *)\n",
|
|
"let m_b12 = Mat.init_cols n_ao n_ao (fun i j -> \n",
|
|
"integral_general (fun a b e f i j ->\n",
|
|
" ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}\n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
"\n",
|
|
"(* Calcul de toutes les intégrales A_12 -> Matrice *)\n",
|
|
"let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->\n",
|
|
"integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j} *. m_C.{e,i} *. m_C.{f,j} \n",
|
|
" -. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} ) \n",
|
|
" *. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )\n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
"\n",
|
|
"(* Calcul de tous les alpha -> Matrice *)\n",
|
|
"let m_alpha = Mat.init_cols n_mo n_mo ( fun i j -> if i= j \n",
|
|
" then 0. \n",
|
|
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))))\n",
|
|
"\n",
|
|
"(*\n",
|
|
"let s_D = \n",
|
|
" let v_D = \n",
|
|
" let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
|
|
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
" in Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
|
|
"in Vec.sum v_D \n",
|
|
"*)\n",
|
|
"\n",
|
|
"(* Calcul de D *)\n",
|
|
"let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
|
|
"integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
|
|
" ) i j\n",
|
|
" )\n",
|
|
" \n",
|
|
"\n",
|
|
"let v_D = Vec.init n_mo ( fun i -> m_D.{i,1} )\n",
|
|
"\n",
|
|
"let s_D = Vec.sum v_D\n",
|
|
"\n",
|
|
"\n",
|
|
"(* \n",
|
|
"(* Calcul de D *)\n",
|
|
"let v_D = Vec.init n_ao (fun i ->\n",
|
|
"integral_general (fun a b e f i _ -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,i} \n",
|
|
" )\n",
|
|
" )\n",
|
|
"let s_D = Vec.sum v_D\n",
|
|
"*)\n",
|
|
"\n",
|
|
"(* Rotation des orbitales *)\n",
|
|
"\n",
|
|
"(* Extraction du plus grand angle alpha*)\n",
|
|
"let ij_max_alpha = Mat.as_vec m_alpha\n",
|
|
"|> iamax;;\n",
|
|
"\n",
|
|
"(* Indices du plus grand angle alpha *)\n",
|
|
"let indice_i = (ij_max_alpha - 1) mod n_ao + 1;;\n",
|
|
"let indice_j = (ij_max_alpha - 1) / n_mo + 1;;\n",
|
|
"\n",
|
|
"(* Valeur du plus grand angle alpha*)\n",
|
|
"let alpha = m_alpha.{indice_i,indice_j};;\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"let m_Phi = m_C\n",
|
|
"\n",
|
|
"(* Matrice de rotation 2 par 2 *)\n",
|
|
"let m_R = Mat.init_cols 2 2 (fun i j -> if i=j then cos alpha\n",
|
|
" else if i>j then sin alpha \n",
|
|
" else -. sin alpha );;\n",
|
|
"\n",
|
|
"(* 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)\n",
|
|
"pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n",
|
|
"let p = indice_i;;\n",
|
|
"let q = indice_j;;\n",
|
|
"let m_Ksi = Mat.init_cols n_ao 2 (fun i j -> if j=1 then m_Phi.{i,p} else m_Phi.{i,q} );;\n",
|
|
"\n",
|
|
"(* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle on\n",
|
|
"obtient les deux orbitales que l'on va réinjecter dans la matrice Phi*)\n",
|
|
"let m_Ksi_tilde = gemm m_Ksi m_R;;\n",
|
|
"\n",
|
|
"\n",
|
|
"(* Réinjection*)\n",
|
|
"\n",
|
|
"\n",
|
|
"(* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur \n",
|
|
"les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice\n",
|
|
"des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec \n",
|
|
"celle contenant i~ et j~ *)\n",
|
|
"\n",
|
|
"(* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n",
|
|
"let m_Psi_tilde = Mat.init_cols n_ao n_ao (fun i j -> if j=q then m_Ksi_tilde.{i,2}\n",
|
|
" else if j=p then m_Ksi_tilde.{i,1}\n",
|
|
" else 0.);;\n",
|
|
" (* p q verif*) \n",
|
|
"(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
|
|
"let m_Psi = Mat.init_cols n_ao n_ao (fun i j -> if j=q then m_Ksi.{i,2}\n",
|
|
" else if j=p then m_Ksi.{i,1}\n",
|
|
" else 0.);;\n",
|
|
" \n",
|
|
"(* 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\n",
|
|
"par la matrice *)\n",
|
|
"let m_interm = Mat.sub m_Phi m_Psi;;\n",
|
|
"\n",
|
|
"(* Construction de la nouvelle matrice d'OMs phi~ *)\n",
|
|
"let m_Phi_tilde = Mat.add m_Psi_tilde m_interm;; \n",
|
|
"\n",
|
|
"*)"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "OCaml 4.10.0",
|
|
"language": "OCaml",
|
|
"name": "ocaml-jupyter"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": "text/x-ocaml",
|
|
"file_extension": ".ml",
|
|
"mimetype": "text/x-ocaml",
|
|
"name": "OCaml",
|
|
"nbconverter_exporter": null,
|
|
"pygments_lexer": "OCaml",
|
|
"version": "4.10.0"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|