diff --git a/Work.ipynb b/Work.ipynb index 877a68c..72597c3 100644 --- a/Work.ipynb +++ b/Work.ipynb @@ -203,24 +203,24 @@ "metadata": {}, "outputs": [], "source": [ - "(* Fonction création chaîne de n H *)\n", + "(* Fonction création chaîne linéaire de n H *)\n", "\n", - "let xyz n = \n", + "let xyz d n = \n", " let accu = \"\"\n", " in\n", - " let rec toto accu n =\n", + " let rec toto accu d n =\n", " let accu = \n", " if n=0 \n", " then accu ^ \"\"\n", " else\n", " match n with \n", - " | 1 -> \"H\" ^ \" \" ^ string_of_float(1.8 *. float_of_int(n)) ^ \" 0. 0.\\n\" ^ accu\n", - " | x -> toto (\"H\" ^ \" \" ^ string_of_float(1.8 *. float_of_int(n)) ^ \" 0. 0.\\n\" ^ accu ) (n-1)\n", + " | 1 -> \"H 0. 0. 0.\\n\" ^ accu\n", + " | x -> toto (\"H\" ^ \" \" ^ string_of_float( d *. float_of_int(n-1)) ^ \" 0. 0.\\n\" ^ accu ) d (n-1)\n", " in\n", " accu\n", - " in string_of_int(n) ^ \"\\nH\" ^ string_of_int(n) ^ \"\\n\" ^ toto accu n;;\n", + " in string_of_int(n) ^ \"\\nH\" ^ string_of_int(n) ^ \"\\n\" ^ toto accu d n;;\n", " \n", - "let xyz_string = xyz 4;;" + "let xyz_string = xyz 1.8 4;;" ] }, { @@ -1054,7 +1054,15 @@ "\n", "$A^r_{12}=A^x_{12} + A^y_{12} + A^z_{12}$\n", "\n", - "$B^r_{12}=B^x_{12} + B^y_{12} + B^z_{12}$" + "$B^r_{12}=B^x_{12} + B^y_{12} + B^z_{12}$\n", + "\n", + "Et le critère à minimiser est :\n", + "\n", + "$D= \\sum_i < \\phi_i | r | \\phi_i >$\n", + "\n", + "Avec \n", + "\n", + "$< \\phi_i | r | \\phi_i > = < \\phi_i | x | \\phi_i > + < \\phi_i | y | \\phi_i > + < \\phi_i | z | \\phi_i >$\n" ] }, { @@ -1247,7 +1255,8 @@ "\n", "\n", "Avec : \n", - " $A_{ij} = [ \\phi_i \\phi_j | \\phi_i \\phi_j] - 1/4 [\\phi_i^2 - \\phi_j^2 | \\phi_i^2 - \\phi_j^2] $\n", + "\n", + "$A_{ij} = [ \\phi_i \\phi_j | \\phi_i \\phi_j] - 1/4 [\\phi_i^2 - \\phi_j^2 | \\phi_i^2 - \\phi_j^2] $\n", " \n", "Où :\n", "\n", @@ -1627,8 +1636,9 @@ " then 0. \n", " else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))));;\n", "\n", - "f_alpha m_C;; \n", - " " + "(*********************)\n", + "\n", + "f_alpha m_C;; " ] }, { @@ -1644,11 +1654,12 @@ " 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.init n_ao ( fun i -> m_D.{i,i} )\n", " in Vec.sum v_D ;;\n", "\n", + "(******************)\n", "\n", - "s_D m_C;;\n" + "s_D m_C;;" ] }, { @@ -1669,7 +1680,7 @@ " |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef \n", " in\n", " let phi_z_phi =\n", - " Multipole.matrix_y multipoles \n", + " Multipole.matrix_z multipoles \n", " |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef\n", " in \n", "\n", @@ -1689,10 +1700,54 @@ " Mat.init_cols n_ao n_ao ( fun i j ->\n", " 0.25 *. asin(m_b12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) )));;\n", "\n", + "(****************************)\n", "\n", "f_alpha_boys multipoles;;\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(* Fonction de calcul de D Boys *)\n", + "\n", + "let d_boys multipoles = \n", + "\n", + " let phi_x2_phi =\n", + " Multipole.matrix_x2 multipoles \n", + " |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef \n", + " in\n", + " \n", + " (*Util.debug_matrix \"phi_x2_phi\" phi_x2_phi;*)\n", + " \n", + " let phi_y2_phi =\n", + " Multipole.matrix_y2 multipoles \n", + " |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef \n", + " in\n", + " \n", + " (*Util.debug_matrix \"phi_y2_phi\" phi_y2_phi;*)\n", + " \n", + " let phi_z2_phi =\n", + " Multipole.matrix_z2 multipoles \n", + " |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef\n", + "\n", + " in\n", + " \n", + " (*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", + "\n", + "in\n", + "Vec.sum v_D_boys;;\n", + "\n", + "(*************************)\n", + "\n", + "d_boys multipoles;;" + ] + }, { "cell_type": "code", "execution_count": null, @@ -1717,6 +1772,8 @@ "in \n", "alpha methode;;\n", "\n", + "(*************************)\n", + "\n", "calcul_alpha_2 \"ER\" ;;\n", "calcul_alpha_2 \"Boys\" ;;\n", "let methode = \"ER\";;\n", @@ -1732,7 +1789,7 @@ "(* Test méthode de calcul du critère de convergence *)\n", "let calcul_cc methode = \n", "\n", - "let b_boys = 2.\n", + "let b_boys = d_boys multipoles\n", "in\n", "let d_er = s_D m_C\n", "in\n", @@ -1747,7 +1804,9 @@ "in \n", "cc methode;;\n", "\n", - "calcul_cc \"ER\";;\n" + "(*************************)\n", + "\n", + "calcul_cc \"boys\";;" ] }, { @@ -1821,6 +1880,8 @@ " then {alpha_max; indice_ii; indice_jj}\n", " else new_m_alpha alpha_m (n_rec_alpha-1);;\n", "\n", + "(*************************)\n", + "\n", "(*\n", "let m_alpha = f_alpha m_C;;\n", "new_m_alpha m_alpha 3;;\n", @@ -1842,6 +1903,8 @@ " then sin alpha \n", " else -. sin alpha );;\n", "\n", + "(*************************)\n", + "\n", "(*\n", "let alpha = alphaij.alpha_max;; (* Fonction -> constante *) \n", "f_R alpha;;\n", @@ -1876,6 +1939,8 @@ "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", + "\n", "(*\n", "let m_Ksi = f_Ksi indice_i indice_j m_C;; (* Fonction -> constante *)\n", "*)" @@ -1890,6 +1955,9 @@ "(* 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", + "(*************************)\n", + "\n", "(*\n", "let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi;; (* Fonction -> constante *)\n", "*)" @@ -1922,6 +1990,9 @@ " then m_Ksi.{i,2}\n", " else 0.)\n", ";;\n", + "\n", + "(*************************)\n", + "\n", "(*\n", "let m_Psi = f_Psi m_Ksi indice_i indice_j;; (* Fonction -> constante *)\n", "\n", @@ -1938,6 +2009,9 @@ "(* 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", + "(*************************)\n", + "\n", "(*\n", "let m_interm = f_interm m_C m_Psi;; (* Fonction -> constante *)\n", "\n",