Work : correction et mise à jour des fonctions, matrice partielle ok

This commit is contained in:
Yann Damour 2020-05-15 12:33:33 +02:00
parent 5eff4a834c
commit c478f9ed4e

View File

@ -161,13 +161,15 @@
},
"outputs": [],
"source": [
"(*\n",
"let xyz_string = \n",
"\"3\n",
"Water\n",
"O 0. 0. 0.\n",
"H -0.756950272703377558 0. -0.585882234512562827\n",
"H 0.756950272703377558 0. -0.585882234512562827\n",
"\""
"\"\n",
"*)"
]
},
{
@ -177,24 +179,6 @@
"### H$_4$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"let xyz_string = \n",
"\"4\n",
"H4\n",
"H 0. 0. 0.\n",
"H 1.8 0. 0.\n",
"H 3.6 0. 0.\n",
"H 5.4 0. 0.\n",
"\"\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
@ -202,7 +186,7 @@
"outputs": [],
"source": [
"(* Fonction création chaîne linéaire de n H *)\n",
"(*\n",
"\n",
"let xyz d n = \n",
" let accu = \"\"\n",
" in\n",
@ -218,8 +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 4;;\n"
]
},
{
@ -250,7 +233,7 @@
},
"outputs": [],
"source": [
"(*\n",
"\n",
"let basis_string = \n",
"\"\n",
"HYDROGEN\n",
@ -263,10 +246,10 @@
"6 0.1001124280E+00 0.1303340841E+00\n",
"\n",
"\"\n",
"*)\n",
"\n",
"\n",
"\n",
"(*\n",
"let basis_string = \n",
"\"\n",
"HYDROGEN\n",
@ -312,7 +295,8 @@
"1 2.753000E-01 1.000000E+00\n",
"D 1\n",
"1 1.185000E+00 1.0000000\n",
"\""
"\"\n",
"*)"
]
},
{
@ -1497,7 +1481,7 @@
"type alphaij = {\n",
" alpha_max : float;\n",
" indice_ii : int;\n",
" indice_jj : int;};;"
" indice_jj : int;};;\n"
]
},
{
@ -1506,6 +1490,8 @@
"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",
@ -1523,7 +1509,39 @@
" in sum v\n",
") (Util.array_range 1 n_ao)\n",
"|> sum \n",
"\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",
"*)"
]
},
{
@ -1532,6 +1550,99 @@
"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",
" \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;;"
]
},
{
"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",
@ -1567,43 +1678,8 @@
"\n",
"(*********************)\n",
"\n",
"f_alpha m_C 1.e-5;; "
]
},
{
"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",
" 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;; "
"f_alpha m_C 1.e-5;; \n",
"*)"
]
},
{
@ -1738,7 +1814,6 @@
" in\n",
" let d_boys = d_boys m_C\n",
" in\n",
" (*Printf.printf \"%f\\n\" d_boys;*)\n",
" let alpha_er = f_alpha m_C\n",
" in\n",
" let d_er = s_D m_C\n",
@ -1795,8 +1870,9 @@
"source": [
"(* 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 n_mo = Mat.dim1 m_alpha\n",
"let rec new_m_alpha m_alpha m_C n_rec_alpha=\n",
"\n",
" let n_mo = Mat.dim2 m_C\n",
" in\n",
" let alpha_m =\n",
" \n",
@ -1858,13 +1934,14 @@
" \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",
" else new_m_alpha alpha_m m_C (n_rec_alpha-1);;\n",
"\n",
"(*************************)\n",
"\n",
"\n",
"let m_alpha = f_alpha m_C;;\n",
"new_m_alpha m_alpha 3;;\n"
"(*\n",
"let m_alpha = f_alpha occ\n",
"let alphaij = new_m_alpha m_alpha occ 3;;\n",
"alphaij.alpha_max;;\n",
"*)"
]
},
{
@ -1916,13 +1993,16 @@
"source": [
"(* 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",
"let f_Ksi indice_i indice_j m_C =\n",
" let n_ao = Mat.dim1 m_C\n",
"in\n",
"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",
"(*************************)\n",
"\n",
"(*\n",
"let m_Ksi = f_Ksi indice_i indice_j m_C;; (* Fonction -> constante *)\n",
"*)"
"*)\n"
]
},
{
@ -1933,7 +2013,7 @@
"source": [
"(* 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",
"let f_Ksi_tilde m_R m_Ksi m_C = gemm m_Ksi m_R;;\n",
"\n",
"(*************************)\n",
"\n",
@ -1954,7 +2034,12 @@
"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_mo (fun i j -> \n",
"let f_Psi_tilde m_Ksi_tilde indice_i indice_j m_C=\n",
" let n_mo = Mat.dim2 m_C\n",
" in\n",
" let n_ao = Mat.dim1 m_C\n",
" in\n",
"Mat.init_cols n_ao n_mo (fun i j -> \n",
" if j=indice_i \n",
" then m_Ksi_tilde.{i,1}\n",
" else if j=indice_j then m_Ksi_tilde.{i,2}\n",
@ -1962,7 +2047,13 @@
";;\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_mo (fun i j -> \n",
"let f_Psi m_Ksi indice_i indice_j m_C = \n",
" let n_mo = Mat.dim2 m_C\n",
" in\n",
" let n_ao = Mat.dim1 m_C\n",
" in\n",
" Printf.printf \"%i %i\\n\" n_mo n_ao;\n",
"Mat.init_cols n_ao n_mo (fun i j -> \n",
" if j=indice_i \n",
" then m_Ksi.{i,1}\n",
" else if j=indice_j \n",
@ -2007,7 +2098,7 @@
"outputs": [],
"source": [
"let toto = [4];;\n",
"let occ_m_C m_C toto= Mat.init_cols n_ao 3 (fun i j ->\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",
@ -2067,7 +2158,7 @@
" (* 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 n_rec_alpha (* Fonction -> constante *)\n",
" let alphaij = new_m_alpha m_alpha m_C n_rec_alpha (* Fonction -> constante *)\n",
" in\n",
"\n",
" (* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)\n",
@ -2095,19 +2186,19 @@
" (*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",
" let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi m_C (* 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",
" let m_Psi = f_Psi m_Ksi indice_i indice_j m_C(* 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",
" let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j m_C (* Fonction -> constante *)\n",
" in\n",
"\n",
" (*Util.debug_matrix \"m_Psi_tilde\" m_Psi_tilde;*)\n",
@ -2143,10 +2234,12 @@
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"(* Calcul *)\n",
"(* Fonction / Matrice des coef / Méthode(\"Boys\" ou \"ER\") / 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 \"boys\" 1. 20 0. 10e-7;;"
"final_m_C m_C \"ER\" 1. 1 0. 10e-7;;\n",
"*)"
]
},
{
@ -2155,19 +2248,7 @@
"metadata": {},
"outputs": [],
"source": [
"let new_m_C = final_m_C occ \"ER\" 1. 50 0. 10e-15;;\n",
"let m_alpha = f_alpha new_m_C;;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m_alpha_d \"Boys\" new_m_C;;\n",
"d_boys m_C;;\n",
"d_boys new_m_C;;"
"let new_m_C = final_m_C occ \"ER\" 1. 20 0. 10e-15;;\n"
]
},
{