From c478f9ed4eaf50c64b12af8273df857fd1a29c59 Mon Sep 17 00:00:00 2001 From: yann Date: Fri, 15 May 2020 12:33:33 +0200 Subject: [PATCH] =?UTF-8?q?Work=20:=20correction=20et=20mise=20=C3=A0=20jo?= =?UTF-8?q?ur=20des=20fonctions,=20matrice=20partielle=20ok?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Work.ipynb | 273 ++++++++++++++++++++++++++++++++++------------------- 1 file changed, 177 insertions(+), 96 deletions(-) diff --git a/Work.ipynb b/Work.ipynb index b49c4b0..4929617 100644 --- a/Work.ipynb +++ b/Work.ipynb @@ -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) = = \\sum_p

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) = = \\sum_q 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) = = \\sum_r 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" ] }, {