Work : update test m_d, test match lo/deloc= echec

This commit is contained in:
Yann Damour 2020-05-15 16:31:00 +02:00
parent c478f9ed4e
commit a68d34aca0

View File

@ -202,7 +202,7 @@
" accu\n",
" in string_of_int(n) ^ \"\\nH\" ^ string_of_int(n) ^ \"\\n\" ^ toto accu d n;;\n",
" \n",
"let xyz_string = xyz 1.8 4;;\n"
"let xyz_string = xyz 1.8 6;;\n"
]
},
{
@ -244,7 +244,6 @@
"4 0.6259552659E+00 0.3705627997E+00\n",
"5 0.2430767471E+00 0.4164915298E+00\n",
"6 0.1001124280E+00 0.1303340841E+00\n",
"\n",
"\"\n",
"\n",
"\n",
@ -1484,66 +1483,6 @@
" indice_jj : int;};;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\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",
"(* 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,
@ -1633,7 +1572,174 @@
"\n",
"(*********************)\n",
"\n",
"f_alpha m_C;;"
"f_alpha m_C;;\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\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",
"(* 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": [
"(* 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 =\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",
" let m_d = 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",
" 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",
" m_d.{i,j} <- m_d.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j}\n",
" ) (Util.array_range 1 n_mo) \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",
" ),Mat.init_cols n_mo n_mo ( fun i j -> m_d.{i,j}));;\n",
"\n",
"(*********************)\n",
"\n",
"f_alpha m_C;;\n",
"\n",
"let m_alpha , m_d = f_alpha m_C;;\n",
"\n",
"\n",
"let s_D m_C = \n",
" let v_D = \n",
" let n_mo = Mat.dim2 m_C\n",
" in Vec.init n_mo ( fun i -> m_d.{i,i} )\n",
"in Vec.sum v_D ;;\n",
"\n",
"\n",
"s_D m_C;;\n"
]
},
{
@ -1688,6 +1794,7 @@
"metadata": {},
"outputs": [],
"source": [
"\n",
"(* Calcul de D -> critère à maximiser dans ER*)\n",
"let s_D m_C = \n",
" let v_D = \n",
@ -1701,8 +1808,12 @@
" in Vec.sum v_D ;;\n",
"\n",
"(******************)\n",
"\n",
"s_D m_C;;"
"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",
"toto *. 6.;;\n"
]
},
{
@ -1747,8 +1858,9 @@
" 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_boys m_C;;"
"(*\n",
"f_alpha_boys m_C;;\n",
"*)"
]
},
{
@ -1937,11 +2049,10 @@
" else new_m_alpha alpha_m m_C (n_rec_alpha-1);;\n",
"\n",
"(*************************)\n",
"(*\n",
"let m_alpha = f_alpha occ\n",
"let alphaij = new_m_alpha m_alpha occ 3;;\n",
"alphaij.alpha_max;;\n",
"*)"
"\n",
"let m_alpha = f_alpha m_C\n",
"let alphaij = new_m_alpha m_alpha m_C 3;;\n",
"alphaij.alpha_max;;\n"
]
},
{
@ -1950,21 +2061,43 @@
"metadata": {},
"outputs": [],
"source": [
"let f_id loc_deloc =\n",
" let id_loc = Mat.identity 2\n",
" in \n",
" let id_deloc = Mat.init_cols 2 2 (fun i j -> \n",
" if i = j \n",
" then -.1.\n",
" else 0. )\n",
" in\n",
" let toto loc_deloc = \n",
" match loc_deloc with \n",
" | \"loc\" -> id_loc\n",
" | \"deloc\" -> id_deloc\n",
" | _ -> invalid_arg \"Unknown method, please enter loc or deloc\"\n",
" \n",
"in toto loc_deloc;;\n",
"\n",
"f_id \"deloc\";;\n",
"\n",
"(* Matrice de rotation 2 par 2 *)\n",
"let f_R alpha =\n",
" Mat.init_cols 2 2 (fun i j -> \n",
" if i=j \n",
" then cos alpha\n",
" else if i>j \n",
" then sin alpha \n",
" else -. sin alpha );;\n",
"let f_R alpha loc_deloc =\n",
" let m_id = f_id loc_deloc\n",
" in\n",
" let m_R = \n",
" Mat.init_cols 2 2 (fun i j -> \n",
" if i=j \n",
" then cos alpha\n",
" else if i>j \n",
" then sin alpha \n",
" else -. sin alpha )\n",
" in gemm m_id m_R;;\n",
"\n",
"(*************************)\n",
"\n",
"(*\n",
"\n",
"let alpha = alphaij.alpha_max;; (* Fonction -> constante *) \n",
"f_R alpha;;\n",
"*)"
"f_R alpha \"loc\";;\n",
"f_R alpha \"deloc\";;\n"
]
},
{
@ -2019,7 +2152,7 @@
"\n",
"(*\n",
"let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi;; (* Fonction -> constante *)\n",
"*)"
"*)\n"
]
},
{
@ -2097,13 +2230,60 @@
"metadata": {},
"outputs": [],
"source": [
"let toto = [4];;\n",
"(* Test localisation matrice rectangulaire et partielle *)\n",
"(*let toto = [4];;\n",
"let occ_m_C m_C toto= Mat.init_cols 4 3 (fun i j ->\n",
" if not (List.mem j toto) \n",
" then m_C.{i,j}\n",
" else 0.);;\n",
" \n",
"let occ = occ_m_C m_C toto;;"
"let occ = occ_m_C m_C toto;;\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",
"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",
" \n",
" in\n",
" List.filter (fun x -> x > 0) vec_list_int;;\n",
"\n",
"(* 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\n",
"dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)\n",
"let split_mat mat list =\n",
" let vec_of_mat = Mat.to_col_vecs mat\n",
" in\n",
" let f a = vec_of_mat.(a-1)\n",
" in\n",
" let vec_list_1 = List.map f list\n",
" in\n",
" let list_2 = miss_elem mat list\n",
" in\n",
" let vec_list_2 = List.map f list_2\n",
"in (Mat.of_col_vecs_list vec_list_1,Mat.of_col_vecs_list vec_list_2);;\n",
"\n",
"m_C;;\n",
"\n",
"let list_om = [1;2]\n",
"\n",
"let m_occ , m_vir = split_mat m_C list_om;;"
]
},
{
@ -2122,7 +2302,7 @@
"(* Localisation de Edminstion ou de Boys *)\n",
"\n",
"(* Calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)\n",
"let rec final_m_C m_C methode epsilon n prev_critere_D cc=\n",
"let rec final_m_C m_C methode loc_deloc epsilon n prev_critere_D cc=\n",
"\n",
" Printf.printf \"%i\\n%!\" n;\n",
"\n",
@ -2133,7 +2313,7 @@
" 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",
" let new_m_C m_C methode loc_deloc =\n",
" \n",
" (* Fonction de pattern matching en fonction de la méthode *)\n",
" let alphad = m_alpha_d methode m_C \n",
@ -2165,6 +2345,8 @@
" let alpha = epsilon *. alphaij.alpha_max (* Fonction -> constante *)\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 (* Fonction -> constante *)\n",
" in\n",
@ -2172,12 +2354,12 @@
" 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",
" \n",
" (* Matrice de rotation *)\n",
" let m_R = f_R alpha loc_deloc (* Fonction -> constante *)\n",
" in\n",
"\n",
" (*Util.debug_matrix \"m_R\" m_R;*)\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",
@ -2212,7 +2394,7 @@
" (* Matrice après rotation *)\n",
" ( Mat.add m_Psi_tilde m_interm, critere_D)\n",
" in\n",
" let m_new_m_C , critere_D = new_m_C m_C methode (* Fonction -> constante *)\n",
" let m_new_m_C , critere_D = new_m_C m_C methode loc_deloc(* Fonction -> constante *)\n",
" in\n",
" let diff = prev_critere_D -. critere_D +. 1.\n",
" \n",
@ -2225,7 +2407,7 @@
" then m_new_m_C\n",
" else\n",
"\n",
"final_m_C m_new_m_C methode epsilon (n-1) critere_D cc;;"
"final_m_C m_new_m_C methode loc_deloc epsilon (n-1) critere_D cc;;"
]
},
{
@ -2234,12 +2416,10 @@
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"(* Calcul *)\n",
"(* Fonction / Matrice des coef / Méthode(\"Boys\" ou \"ER\") / Pas(=< 1.) / Nombre d'itérations max / \n",
"(* Fonction / Matrice des coef / Méthode(\"Boys\" ou \"ER\") / Localisation ou non (\"loc\" ou \"deloc\"/ Pas(<=1.) / Nombre d'itérations max / \n",
"0. (valeur de D pour initier la boucle) / critère de convergence sur D*)\n",
"final_m_C m_C \"ER\" 1. 1 0. 10e-7;;\n",
"*)"
"let new_m = final_m_C m_C \"ER\" \"loc\" 1. 5 0. 10e-7;;\n"
]
},
{
@ -2248,7 +2428,7 @@
"metadata": {},
"outputs": [],
"source": [
"let new_m_C = final_m_C occ \"ER\" 1. 20 0. 10e-15;;\n"
"final_m_C new_m \"ER\" \"deloc\" 1. 5 0. 10e-7;;\n"
]
},
{