From 514ea607e2a020afbc1c48dbbafc8e320124847c Mon Sep 17 00:00:00 2001 From: yann Date: Thu, 7 May 2020 18:00:34 +0200 Subject: [PATCH] =?UTF-8?q?Work:=20r=C3=A9=C3=A9criture=20boucle=20ER=20av?= =?UTF-8?q?ec=20fct=20def=20ext=20et=20propre?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Work.ipynb | 845 ++++++++++++++++++++++++++++------------------------- 1 file changed, 455 insertions(+), 390 deletions(-) diff --git a/Work.ipynb b/Work.ipynb index 4063908..5993aea 100644 --- a/Work.ipynb +++ b/Work.ipynb @@ -186,10 +186,12 @@ "outputs": [], "source": [ "let xyz_string = \n", - "\"2\n", - "H2\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", "\"" ] }, @@ -221,18 +223,22 @@ }, "outputs": [], "source": [ - "let basis_string = \"\n", - "HYDROGEN\n", - "S 4\n", - "1 1.301000E+01 1.968500E-02\n", - "2 1.962000E+00 1.379770E-01\n", - "3 4.446000E-01 4.781480E-01\n", - "4 1.220000E-01 5.012400E-01\n", - "S 1\n", - "1 1.220000E-01 1.000000E+00\n", - "P 1\n", - "1 7.270000E-01 1.0000000\n", + "let basis_string = \n", "\"\n", + "HYDROGEN\n", + "S 6\n", + "1 0.3552322122E+02 0.9163596281E-02\n", + "2 0.6513143725E+01 0.4936149294E-01\n", + "3 0.1822142904E+01 0.1685383049E+00\n", + "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", + "\n", + "\n", "(*\n", "let basis_string = \"\n", "HYDROGEN\n", @@ -1302,56 +1308,6 @@ "## Calculs : Localisation de Edminston" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "(* Test algo alpha < pi/. 2. *)\n", - "(*)\n", - "let m_alpha = f_alpha m_Phi\n", - "in\n", - "let prev_alpha = m_alpha\n", - "in\n", - "let rec m_new_alpha m_alpha n = \n", - " \n", - " Printf.printf \"%d\\n%!\" n;\n", - " if n == 0 then\n", - " prev_alpha\n", - " else\n", - " \n", - " \n", - " let m_alpha =\n", - "\n", - " Mat.init_cols n_ao n_ao (fun i j -> if m_alpha.{i,j} > (3.14159265359 /. 2.) \n", - " then (m_alpha.{i,j} -. (3.14159265359 /. 2.))\n", - " else m_alpha.{i,j})\n", - "\n", - " in\n", - " (* Extraction du plus grand angle alpha*)\n", - " let ij_max_alpha = Mat.as_vec m_alpha\n", - " |> iamax\n", - " in\n", - " (* Indices du plus grand angle alpha *)\n", - " let indice_i = ij_max_alpha mod n_ao\n", - " in\n", - " let indice_j = ij_max_alpha / n_ao + 1\n", - " in\n", - " (* Valeur du plus grand angle alpha*)\n", - " let alpha = m_alpha.{indice_i,indice_j}\n", - " in\n", - "\n", - " if alpha < (3.14159265359 /. 2.)\n", - " then\n", - " m_alpha\n", - " else\n", - " let m_alpha = prev_alpha\n", - " in\n", - " m_new_alpha m_alpha (n-1);;\n", - "*)" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1443,10 +1399,9 @@ " )\n", "\n", "(* Calcul de tous les alpha -> Matrice *)\n", - "let m_alpha = Mat.init_cols n_mo n_mo ( fun i j ->\n", - " asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\n", - " )\n", - " )\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", @@ -1552,8 +1507,22 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Réécriture pour itérerer ER\n", - "-> Il fonctionne apparemment mais je suis sceptique pour certains résultats, il faut réécrire toutes les fonctions en fonction de leur arguments pour le refaire mieux et éviter les problèmes rencontrés sur les variables locales et globales" + "## Réécriture pour itérerer ER" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Réécriture des fonctions avec leurs arguments \n", + "-> pour refaire mieux la boucle en évitant des problèmes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Fonctions" ] }, { @@ -1562,19 +1531,27 @@ "metadata": {}, "outputs": [], "source": [ - "(*\n", - "\n", - "(* Fonction pour itérer , cc critère de convergence*) \n", - "(* Localisation de Edminstion *)\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", - "(*Calcul des intégrales*) \n", + " \n", + "type alphaij = {\n", + " alpha_max : float;\n", + " indice_ii : int;\n", + " indice_jj : int;};;" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(*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", @@ -1590,177 +1567,259 @@ " ) (Util.array_range 1 n_ao)\n", " in sum v\n", ") (Util.array_range 1 n_ao)\n", - "|> sum \n", + "|> sum " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* 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", - "let m_Phi = None;;\n", - "type iteration = \n", - "{ coef : Mat.t ;\n", - " loc : float ;\n", - " }\n", - "\n", - "(*let prev_iter = {coef = m_C ;loc = 0. }*)\n", - "\n", - "(* Calcul de la nouvelle matrice Phi après rotations *)\n", - "let rec new_Phi prev_iter n = \n", - " Printf.printf \"%d\\n%!\" n;\n", - " \n", - " Util.debug_matrix \"prev_iter\" prev_iter.coef;\n", - " \n", - " if n == 0 then\n", - " prev_iter\n", - " else\n", - " let m_Phi = prev_iter.coef\n", - " in\n", - " \n", - " Util.debug_matrix \"m_coef\" prev_iter.coef;\n", - " (*Util.debug_matrix \"m_Phi\" m_Phi;\n", - " *)\n", - " (* Calcul de tous les alpha -> Matrice *)\n", - " let f_alpha m_Phi =\n", - "\n", - " (* Calcul de toutes les intégrales B_12 -> Matrice *)\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_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} ) *. m_Phi.{e,i} *. m_Phi.{f,j}\n", - " ) i j\n", - " )\n", - " in\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_Phi.{a,i} *. m_Phi.{b,j} *. m_Phi.{e,i} *. m_Phi.{f,j} \n", - " -. 0.25 *. ( m_Phi.{e,i} *. m_Phi.{f,i} -. m_Phi.{e,j} *. m_Phi.{f,j} ) \n", - " *. ( m_Phi.{a,i} *. m_Phi.{b,i} -. m_Phi.{a,j} *. m_Phi.{b,j} )\n", - " ) i j\n", - " )\n", - " in\n", - " Mat.init_cols n_mo n_mo ( fun i j ->\n", - " asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )\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", - " (* Calcul de D *)\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_Phi.{a,i} *. m_Phi.{b,i} *. m_Phi.{e,i} *. m_Phi.{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", - " in \n", - " (* Matrice de rotation 2 par 2 *)\n", - " let indice_i, indice_j, m_R = \n", "\n", - " \n", - " let m_alpha = f_alpha m_Phi\n", - " in\n", - " \n", - " Util.debug_matrix \"m_alpha\" m_alpha;\n", - " \n", - " (* Valeur du plus grand angle alpha*)\n", - " (* Extraction du plus grand angle alpha*)\n", - " let ij_max_alpha = Mat.as_vec m_alpha\n", - " |> iamax\n", - " in\n", - " \n", - " Printf.printf \"%i\\n%!\" ij_max_alpha;\n", - " \n", - " (* Indices du plus grand angle alpha *)\n", - " let indice_i = (ij_max_alpha - 1) mod n_mo + 1\n", - " in\n", - " let indice_j = (ij_max_alpha - 1) / n_mo + 1\n", - "\n", - " in\n", - " \n", - " Printf.printf \"%i\\n%!\" indice_i;\n", - " Printf.printf \"%i\\n%!\" indice_j;\n", - " \n", - " let alpha = \n", - " m_alpha.{indice_i,indice_j}\n", - " in\n", - " \n", - " \n", - " Printf.printf \"%f\\n%!\" alpha;\n", - " \n", - " indice_i, indice_j, \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", - " in\n", - " (*\n", - " Util.debug_matrix \"m_R\" m_R;*)\n", - " Printf.printf \"%d %d\\n%!\" indice_i indice_j;\n", - " \n", - " let p = indice_i\n", - " in\n", - " let q = indice_j\n", - " in\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, m_Ksi_tilde = \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", - " (n par 2) pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n", - " Printf.printf \"%d %d\\n%!\" indice_i 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", - " in\n", - " m_Ksi, gemm m_Ksi m_R\n", - "\n", - " in\n", - " Printf.printf \"%d %d\\n%!\" p q;\n", - " Util.debug_matrix \"m_Ksi\" m_Ksi;\n", - " Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n", - " \n", - " (* Construction de la nouvelle matrice d'OMs phi~ *)\n", - " let m_Phi_tilde m_Ksi m_Ksi_tilde = \n", - "\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_mo (fun i j -> if j=p then m_Ksi_tilde.{i,1}\n", - " else if j=q then m_Ksi_tilde.{i,2}\n", - " else 0.)\n", - " in\n", - " Printf.printf \"%d %d\\n%!\" p q;\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 \n", - " matrice Phi par la matrice *)\n", - " let m_interm = \n", - "\n", - " (* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n", - " let m_Psi = Mat.init_cols n_ao n_mo (fun i j -> if j=p then m_Ksi.{i,1}\n", - " else if j=q then m_Ksi.{i,2}\n", - " else 0.)\n", - "\n", - " in\n", - " (*Util.debug_matrix \"m_Psi\" m_Psi;*)\n", - " Mat.sub m_Phi m_Psi\n", - " in\n", - " (*Util.debug_matrix \"m_interm\" m_interm;*)\n", - " Mat.add m_Psi_tilde m_interm\n", - " in \n", - " let coef = m_Phi_tilde m_Ksi m_Ksi_tilde\n", - " in\n", - " \n", - " Util.debug_matrix \"Nouvelle_matrice_des_coef\" coef;\n", - " \n", - " let loc = s_D \n", - " in\n", - " \n", - " Printf.printf \"%f\\n%!\" loc;\n", - " Printf.printf \"%f\\n%!\" prev_iter.loc;\n", - " \n", - " if abs_float (loc -. prev_iter.loc) < 1.e-5 \n", - " then\n", - " { coef ; loc }\n", - " else\n", - " new_Phi {coef; loc;} (n-1)\n", - "\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", - "new_Phi {coef = m_C ;loc = 0. } 0;;\n", + "f_alpha m_C;; \n", + "*) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* 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", + "(*\n", + "s_D m_C;;\n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "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 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", + " \n", + " (*Printf.printf \"%f\\n%!\" alpha_max;*)\n", + " \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", + "let m_alpha = f_alpha m_C;;\n", + "new_m_alpha m_alpha 3;;\n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* 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", + "\n", + "(*\n", + "let alpha = alphaij.alpha_max;; (* Fonction -> constante *) \n", + "f_R alpha;;\n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(*Uniquement pour pouvoir tester les fonctions après cette cellules*)\n", + "(*\n", + "(* Indice i et j du alpha max après calcul *)\n", + "let indice_i = alphaij.indice_ii;; (* Fonction -> constante *)\n", + " \n", + "let indice_j = alphaij.indice_jj;; (* Fonction -> constante *)\n", + " \n", + " \n", + "let m_R = f_R alpha;; (* Fonction -> constante *) \n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "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", + "\n", + "(*\n", + "let m_Ksi = f_Ksi indice_i indice_j m_C;; (* Fonction -> constante *)\n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "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", + "(*\n", + "let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi;; (* Fonction -> constante *)\n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* 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 -> \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", + " else 0.)\n", + ";;\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 -> \n", + " if j=indice_i \n", + " then m_Ksi.{i,1}\n", + " else if j=indice_j \n", + " then m_Ksi.{i,2}\n", + " else 0.)\n", + ";;\n", + "(*\n", + "let m_Psi = f_Psi m_Ksi indice_i indice_j;; (* Fonction -> constante *)\n", + "\n", + "let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j;; (* Fonction -> constante *)\n", + "*)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* 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", + "(*\n", + "let m_interm = f_interm m_C m_Psi;; (* Fonction -> constante *)\n", + "\n", + "let new_m_C m_C= Mat.add m_Psi_tilde m_interm;;\n", + "\n", + "let m_new_m_C = new_m_C m_C;;\n", "*)" ] }, @@ -1768,8 +1827,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Réécriture des fonctions avec leurs arguments \n", - "-> pour refaire mieux la boucle en évitant des problèmes" + "## Boucle localisation avec fonctions définies à l'extérieur" ] }, { @@ -1777,7 +1835,95 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "(* Localisation de Edminstion *)\n", + "\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);;" + ] }, { "cell_type": "code", @@ -1785,6 +1931,18 @@ "metadata": {}, "outputs": [], "source": [ + "(*Cellule de calcul*)\n", + "final_m_C m_C 10;;" + ] + }, + { + "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", @@ -1838,7 +1996,10 @@ " )\n", "in\n", " Mat.init_cols n_ao n_ao ( fun i j ->\n", - " asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.))))\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", @@ -1859,118 +2020,6 @@ " s_D m_C;;\n", " \n", "(******************************************************************************************************)\n", - " \n", - "(* Fonction d'extraction de la valeur maximum de alpha, avec alpha < pi/2 et extraction des indices de celle-ci *)\n", - "\n", - "(*Fonction d'extraction du plus grand élément d'une matrice*)\n", - "(*\n", - "let element_max m = \n", - "\n", - "Mat.as_vec m\n", - "|> iamax;;\n", - "\n", - "let m_alpha = f_alpha m_C ;;\n", - "\n", - "element_max m_alpha;;\n", - "\n", - "\n", - "let max_element f_alpha m_C = \n", - " let m_alpha = f_alpha m_C\n", - " in \n", - " Mat.as_vec m_alpha\n", - " |> iamax;;\n", - "\n", - "max_element f_alpha m_C;;\n", - "\n", - "(************************)\n", - "\n", - "let numero = element_max m_alpha;;\n", - "(*\n", - "let indice_i numero m n_lignes n_colonnes = \n", - " (numero - 1) mod n_lignes + 1;;\n", - " \n", - "indice_i numero m_alpha 10 10;; *)\n", - "\n", - "let indice_i f_alpha m_C max_element n_lignes = \n", - " let max = max_element f_alpha m_C\n", - " in\n", - " (max - 1) mod n_lignes +1 ;;\n", - "\n", - "indice_i f_alpha m_C max_element n_ao;;\n", - "\n", - "let indice_j f_alpha m_C max_element n_colonnes = \n", - " let max = max_element f_alpha m_C\n", - " in\n", - " (max - 1) / n_colonnes +1 ;;\n", - " \n", - "indice_j f_alpha m_C max_element n_ao;;\n", - "\n", - "*)\n", - "(*\n", - "\n", - "let max_element f_alpha m_C = \n", - " let m_alpha = f_alpha m_C\n", - " in \n", - " Mat.as_vec m_alpha\n", - " |> iamax;;\n", - "\n", - "max_element f_alpha m_C;;\n", - "\n", - "let indice_i f_alpha m_C max_element n_lignes = \n", - " let max = max_element f_alpha m_C\n", - " in\n", - " (max - 1) mod n_lignes +1 ;;\n", - "\n", - "indice_i f_alpha m_C max_element n_ao;;\n", - "\n", - "let indice_j f_alpha m_C max_element n_colonnes = \n", - " let max = max_element f_alpha m_C\n", - " in\n", - " (max - 1) / n_colonnes +1 ;;\n", - " \n", - "indice_j f_alpha m_C max_element n_ao;;\n", - "\n", - "(********************************)\n", - "\n", - "(* Elément max de la matrice f_alpha pour la matrice des coefficient m_C *)\n", - "let max_element2 m_C = \n", - " let m_alpha = f_alpha m_C\n", - " in \n", - " Mat.as_vec m_alpha\n", - " |> iamax;;\n", - "\n", - "max_element2 m_C;;\n", - "\n", - "(* Indice i de l'élément max *)\n", - "let indice_ii m_C n_lignes = \n", - " let max = max_element2 m_C\n", - " in\n", - " (max - 1) mod n_lignes +1 ;;\n", - "\n", - "indice_i f_alpha m_C max_element n_ao;;\n", - "\n", - "(* Indice j de l'élément max *)\n", - "let indice_jj m_C n_colonnes = \n", - " let max = max_element2 m_C\n", - " in\n", - " (max - 1) / n_colonnes +1 ;;\n", - " \n", - "indice_jj m_C n_ao;;\n", - "\n", - "(* Valeur max de alpha *)\n", - "let alpha m_C = \n", - " let m_alpha = f_alpha m_C \n", - " in\n", - " let i = indice_ii m_C n_ao\n", - " in\n", - " let j = indice_jj m_C n_ao\n", - " in\n", - " m_alpha.{i,j};;\n", - "\n", - "alpha m_C;;\n", - "\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", @@ -2095,14 +2144,16 @@ "(* 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" + "m_Phi_tilde m_Psi_tilde m_interm;;\n", + "\n", + "*)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Rassemblement des blocs" + "## Programmes avec fonctions définies à l'intérieurs" ] }, { @@ -2113,6 +2164,10 @@ "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", @@ -2150,6 +2205,15 @@ " \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", @@ -2178,13 +2242,15 @@ " (* 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", + " -. 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", - " asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.))))\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", @@ -2211,13 +2277,13 @@ " 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", + " \n", " Printf.printf \"%i\\n%!\" n_rec_alpha;\n", - " *)\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", + " 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", @@ -2225,9 +2291,9 @@ " then -. m_alpha.{i,j}\n", " else m_alpha.{i,j} )\n", " in \n", - " (*\n", + " \n", " Util.debug_matrix \"alpha_m\" alpha_m;\n", - " *)\n", + " \n", " (* Détermination de l'emplacement du alpha max *)\n", " let max_element3 alpha_m = \n", " Mat.as_vec alpha_m\n", @@ -2281,37 +2347,37 @@ " 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", + " 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", + " \n", " Printf.printf \"%i %i\\n%!\" indice_i indice_j;\n", - " *)\n", + " \n", " let m_R = f_R alpha (* Fonction -> constante *)\n", " in\n", - " \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", + " \n", " Util.debug_matrix \"m_Ksi\" m_Ksi;\n", - " *)\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", + " \n", " Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;\n", - " *)\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", @@ -2332,14 +2398,14 @@ " in\n", " let m_Psi = f_Psi m_Ksi indice_i indice_j (* Fonction -> constante *)\n", " in\n", - " (*\n", + " \n", " Util.debug_matrix \"m_Psi\" m_Psi;\n", - " *)\n", + " \n", " let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j (* Fonction -> constante *)\n", " in\n", - " (*\n", + " \n", " Util.debug_matrix \"m_Psi_tilde\" m_Psi_tilde;\n", - " *)\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", @@ -2350,25 +2416,24 @@ "in\n", "let m_new_m_C = new_m_C m_C (* Fonction -> constante *)\n", "in\n", - "(*\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", + "\n", "final_m_C m_new_m_C (n-1);;\n", "\n", "(*****************************)\n", "\n", - "final_m_C m_C 1000;;\n", - "\n" + "final_m_C m_C 10;;\n", + "\n", + "*)" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "\n", - " \n" + "### test" ] }, {