From 35275f7673f5f96b39affc63b274585bc1fef915 Mon Sep 17 00:00:00 2001 From: yann Date: Wed, 13 May 2020 15:46:48 +0200 Subject: [PATCH] =?UTF-8?q?=20Work=20:=20boucle=20avec=20cc=20et=20modif?= =?UTF-8?q?=20pour=20matrice=20non=20carr=C3=A9es?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Work.ipynb | 83 +++++++++++++++++++++++------------------------------- 1 file changed, 36 insertions(+), 47 deletions(-) diff --git a/Work.ipynb b/Work.ipynb index cb5a12b..88d872a 100644 --- a/Work.ipynb +++ b/Work.ipynb @@ -161,7 +161,7 @@ }, "outputs": [], "source": [ - "(*\n", + "\n", "let xyz_string = \n", "\"3\n", "Water\n", @@ -169,7 +169,7 @@ "H -0.756950272703377558 0. -0.585882234512562827\n", "H 0.756950272703377558 0. -0.585882234512562827\n", "\"\n", - " *) " + " " ] }, { @@ -204,7 +204,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", @@ -220,7 +220,8 @@ " 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;;" + "let xyz_string = xyz 1.8 4;;\n", + "*)" ] }, { @@ -251,7 +252,8 @@ }, "outputs": [], "source": [ - "let basis_string = \n", + "(*let basis_string = \n", + "\n", "\"\n", "HYDROGEN\n", "S 6\n", @@ -264,11 +266,12 @@ "\n", "\"\n", "\n", + "*)\n", "\n", "\n", "\n", - "(*\n", - "let basis_string = \"\n", + "let basis_string = \n", + "\"\n", "HYDROGEN\n", "S 4\n", "1 1.301000E+01 1.968500E-02\n", @@ -312,9 +315,7 @@ "1 2.753000E-01 1.000000E+00\n", "D 1\n", "1 1.185000E+00 1.0000000\n", - "\n", - "\"\n", - "*)" + "\"" ] }, { @@ -1536,7 +1537,7 @@ "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", + " 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", @@ -1544,21 +1545,21 @@ "\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", + " 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_ao n_ao ( fun i j -> \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;; " ] }, { @@ -1570,11 +1571,11 @@ "(* Calcul de D -> critère à maximiser dans ER*)\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", + " 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_ao ( fun i -> m_D.{i,i} )\n", + " in Vec.init n_mo ( fun i -> m_D.{i,i} )\n", " in Vec.sum v_D ;;\n", "\n", "(******************)\n", @@ -1605,19 +1606,19 @@ " in \n", "\n", " let m_b12= \n", - " let b12 g = Mat.init_cols n_ao n_ao ( fun i j ->\n", + " let b12 g = Mat.init_cols n_mo n_mo ( fun i j ->\n", " (g.{i,i} -. g.{j,j}) *. g.{i,j})\n", "\n", " in \n", " Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi))\n", " in \n", " let m_a12 =\n", - " let a12 g = Mat.init_cols n_ao n_ao ( fun i j -> \n", + " let a12 g = Mat.init_cols n_mo n_mo ( fun i j -> \n", " g.{i,j} *. g.{i,j} -. 0.25 *. ((g.{i,i} -. g.{j,j}) *. (g.{i,i} -. g.{j,j})))\n", " in\n", " Mat.add (a12 phi_x_phi) ( Mat.add (a12 phi_y_phi) (a12 phi_z_phi))\n", " in\n", - " Mat.init_cols n_ao n_ao ( fun i j -> \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", @@ -1660,7 +1661,7 @@ " (*Util.debug_matrix \"phi_z2_phi\" phi_z2_phi;*)\n", " \n", " let v_D_boys = \n", - " Vec.init n_ao ( fun i -> phi_x2_phi.{i,i} +. phi_y2_phi.{i,i} +. phi_z2_phi.{i,i})\n", + " Vec.init n_mo ( fun i -> phi_x2_phi.{i,i} +. phi_y2_phi.{i,i} +. phi_z2_phi.{i,i})\n", "\n", "in\n", "Vec.sum v_D_boys;;\n", @@ -1751,7 +1752,7 @@ " \n", " if n_rec_alpha == 0 \n", " then m_alpha \n", - " else Mat.init_cols n_ao n_ao (fun i j -> \n", + " else Mat.init_cols n_mo n_mo (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", @@ -1783,7 +1784,7 @@ " let indice_jj = \n", " let max = max_element3 alpha_m (* Fonction -> constante *)\n", " in\n", - " (max - 1) / n_ao +1\n", + " (max - 1) / n_mo +1\n", " in\n", " \n", " (* Valeur du alpha max*)\n", @@ -1809,10 +1810,9 @@ "\n", "(*************************)\n", "\n", - "(*\n", + "\n", "let m_alpha = f_alpha m_C;;\n", - "new_m_alpha m_alpha 3;;\n", - "*)" + "new_m_alpha m_alpha 3;;\n" ] }, { @@ -1964,9 +1964,9 @@ "(* 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_norme_alpha prev_critere_D =\n", + "let rec final_m_C m_C methode epsilon n prev_critere_D cc=\n", "\n", - " (*Printf.printf \"%i\\n%!\" n;*)\n", + " Printf.printf \"%i\\n%!\" n;\n", "\n", " (*Util.debug_matrix \"m_C\" m_C;*)\n", "\n", @@ -1993,7 +1993,7 @@ " let norme_alpha = norme m_alpha\n", " in\n", " \n", - " (*Printf.printf \"%f\\n%!\" norme_alpha;*)\n", + " Printf.printf \"%f\\n%!\" norme_alpha;\n", " \n", " (*Util.debug_matrix \"m_alpha\" m_alpha;*)\n", "\n", @@ -2051,34 +2051,22 @@ " (*Util.debug_matrix \"m_interm\" m_interm;*)\n", " \n", " (* Matrice après rotation *)\n", - " ( Mat.add m_Psi_tilde m_interm, norme_alpha, critere_D)\n", + " ( Mat.add m_Psi_tilde m_interm, critere_D)\n", " in\n", - " let m_new_m_C , norme_alpha , critere_D = new_m_C m_C methode (* Fonction -> constante *)\n", + " let m_new_m_C , critere_D = new_m_C m_C methode (* Fonction -> constante *)\n", " in\n", " let diff = prev_critere_D -. critere_D\n", " \n", - " (*if epsilon <= 0.05\n", - " then 0.1\n", - " else if norme_alpha <= 0.0001 && epsilon >= 0.2\n", - " then epsilon -. 0.2\n", - " else if norme_alpha <= 0.001 && epsilon >= 0.1\n", - " then epsilon -. 0.1\n", - " else if norme_alpha <= 0.01 && epsilon >= 0.05\n", - " then epsilon -. 0.05\n", - " else epsilon *)\n", - " \n", - " \n", - " \n", "in\n", - "Printf.printf \"%f %f %f %f\\n%!\" epsilon prev_norme_alpha norme_alpha diff;\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. < 0.000001**2.\n", + "if diff**2. < cc**2.\n", " then m_new_m_C\n", " else\n", "\n", - "final_m_C m_new_m_C methode epsilon (n-1) norme_alpha critere_D ;;" + "final_m_C m_new_m_C methode epsilon (n-1) critere_D cc;;" ] }, { @@ -2088,8 +2076,9 @@ "outputs": [], "source": [ "(* Calcul *)\n", - "(* Fonction Matrice des coef Méthode(\"Boys\" ou \"ER\") Pas(=< 1.) Nombre d'itérations *)\n", - "final_m_C m_C \"ER\" 1. 100 10. 0.;;" + "(* 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 \"ER\" 1. 1 0. 10e-7;;" ] }, {