diff --git a/Work.ipynb b/Work.ipynb index 8ca694c..4aa88b1 100644 --- a/Work.ipynb +++ b/Work.ipynb @@ -893,7 +893,7 @@ "- une matrice Psi~ (n par n) avec tous ses éléments nuls sauf les colonnes qui contiennent les OMs i~ et j~\n", "Ainsi, on soustrait à la matrice des OMs Phi la matrice Psi pour supprimer les OMs i et j de celle ci puis on additionne cette nouvelle matrice à la matrice Psi~ pour créer la nouvelle matrice des OMs Phi~ avec i~ et j~.\n", "\n", - "On repart de cette nouvelle matrice Phi~ et on cherche la paire d'orbitale (k,l) ayant le plus grand angle de rotation alpha. Et on procède comme nous l'avons précedemment de manière intérative. Le but étant de maximiser D pour la localisation de Edminston et de minimiser B pour la localisation de Boyls." + "On repart de cette nouvelle matrice Phi~ et on cherche la paire d'orbitale (k,l) ayant le plus grand angle de rotation alpha. Et on procède comme nous l'avons précedemment de manière intérative. Le but étant de maximiser D pour la localisation de Edminston et la localisation de Boyls." ] }, { @@ -1022,7 +1022,7 @@ "Sinon on procède comme Edminson Ruedenberg mais avec les intégrales $A_{12}$ et $B_{12}$ définies comme :\n", "\n", "$A^r_{12} = \\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n", - "$*\\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n", + "$.\\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n", "$- \\frac {1}{4}(\\langle \\phi_1 | \\bar r | \\phi_1 \\rangle $\n", "$- \\langle \\phi_2 | \\bar r | \\phi_2 \\rangle . \\langle \\phi_1 | \\bar r | \\phi_1 \\rangle$\n", "$- \\langle \\phi_2 | \\bar r | \\phi_2 \\rangle)$\n", @@ -1268,6 +1268,194 @@ "$= \\sum_n (\\sum_a \\sum_b \\sum_e \\sum_f c_{an} c_{bn} c_{en} c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Localisation des orbitales moléculaires" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Principe de la rotation par paires d'orbitales.\n", + "\n", + "La rotation des orbitales se fait de la manière suivante :\n", + "\n", + "On extrait de la matrice des coefficients des OMs $\\Phi$, (m_C), les orbitales $i$ et $j$ (vecteurs colonnes) pour former une nouvelle matrice $Ksi$ de dimension (n par 2). \n", + "\n", + "Pour effectuer la rotation des orbitales $i$ et $j$ on utilise la matrice de rotation $R$ (m_R) pour la rotation de 2 orbitales, qui est définie comme la matrice :\n", + "\n", + "$R$ =\n", + "\n", + " ( cos(alpha) -sin(alpha) )\n", + "\n", + " ( sin(alpha) cos(alpha) )\n", + " \n", + "La valeur de $\\alpha$ sera la plus grande de la matrice des alphas *Cf Localisation de Edminston Ruedenberg*\n", + " \n", + "On applique $R$ à $Khi$ : $R Ksi = \\tilde Ksi$ \n", + "\n", + "On obtient $\\tilde Ksi$ la matrice contenant les coefficients des deux nouvelles OMs $\\tilde i$ et $\\tilde j$ obtenues par rotation de $i$ et $j$.\n", + "\n", + "On réinjecte ces deux nouvelles orbitales $\\tilde i$ et $\\tilde j$ à la place des anciennes orbitales $i$ et $j$ dans la matrice des coefficients des OMs $\\Phi$, (m_C), ce qui nous donne une nouvelle matrice des coefficient des OMs $\\tilde \\Phi$, (new_m_C).\n", + "Pour cela on créer des matrices intérmédiaires:\n", + "- une matrice $\\Psi$ (n par n) avec tous ses éléments nuls sauf les colonnes qui contiennent les OMs $i$ et $j$ \n", + "- une matrice $\\tilde \\Psi$ (n par n) avec tous ses éléments nuls sauf les colonnes qui contiennent les OMs $\\tilde i$ et $\\tilde j$\n", + "- une matrice m_interm (n par n) où l'on a soustrait $\\Psi$ à $\\Phi$ pour créer des 0 sur les colonnes des OMs $i$ et $j$ \n", + "\n", + "Ainsi, on soustrait à la matrice des coefficient $\\Phi$ la matrice $\\Psi$ pour supprimer les OMs $i$ et $j$ de celle ci puis on additionne cette nouvelle matrice à la matrice $\\tilde \\Psi$ pour créer la nouvelle matrice des coefficients des OMS $\\tilde \\Phi$, avec $\\tilde i$ et $\\tilde j$.\n", + "\n", + "On repart de cette nouvelle matrice des coefficients $\\tilde \\Phi$ et on cherche la paire d'orbitale $(k,l)$ ayant le plus grand angle de rotation alpha (l'angle $\\alpha$ pour la paire d'orbitale $(i,j)$ étant devenui nul après rotation. Et on procède comme nous l'avons précedemment de manière intérative. Le but étant de maximiser $D$ pour la localisation de Edminston et la localisation de Boyls." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Localisation de Edmiston C. & Ruedenberg K. \n", + "\n", + "(Localisation des orbitales par rotation de paires d'orbitales)\n", + "\n", + "Ref :\n", + "*Edmiston, C., & Ruedenberg, K. (1963). Localized Atomic and Molecular Orbitals.\n", + "Reviews of Modern Physics, 35(3), 457–464. doi:10.1103/revmodphys.35.457*\n", + "\n", + "Méthode :\n", + "\n", + "Le but de cette méthode est de maximiser $D$ :\n", + "\n", + "$D(\\phi)= \\sum_n [\\phi_n^2 | \\phi_n^2 ]$\n", + "\n", + "$= \\sum_n < \\phi^2_n | 1/r_{12} | \\phi^2_n >$\n", + " \n", + " \n", + "\n", + " \n", + "Car selon *J. E. Lennard-Jones and J. A. Pople, Proc. Roy. Soc. (London) A202, 166 (1950)*, on peut générer des orbitales équivalentes et celles ci maximiseront probablement la somme des termes d'auto répulsion orbitalaire $D$.\n", + "\n", + "On va créer des nouvelles orbitales $\\tilde i$ et $\\tilde j$ à partir des orbitales $i$ et $j$ par combinaison linéaire de ces dernières tel que :\n", + "\n", + "$i~ (x) = cos(\\gamma) i(x) + sin(\\gamma) j(x)$\n", + "\n", + "$j~ (x) = -sin(\\gamma) i(x) + cos(\\gamma) j(x)$\n", + "\n", + "On part de la matrices de orbitales moléculaires $\\Phi$ et on cherche la paire d'orbitale $(i,j)$ ayant le plus grand angle de rotation $\\alpha$, avec $\\alpha$ défini comme (si ce dernier est supétieur à $\\frac{\\pi}{2}$ on devra soustraire $\\frac{\\pi}{2}$ aux éléments de la matrices) :\n", + "\n", + "\n", + "$Cos(4 \\alpha)= -A_{12} / (A_{12}^2 + B_{12}^2)^{1/2}$\n", + "\n", + "$\\alpha = (1/4) Acos (-A_{12} / (A_{12}^2 + B_{12}^2)^{1/2})$\n", + "\n", + "$Sin(4 \\alpha)= B_{12} / (A_{12}^2 + B_{12}^2)^{1/2}$\n", + "\n", + "$\\alpha = (1/4) Asin (B_{12} / (A_{12}^2 + B_{12}^2)^{1/2})$\n", + "\n", + "$Tan(4 \\alpha) = -B_{12} / A_{12}$\n", + "\n", + "$\\alpha = (1/4) Atan (-B_{12} / A_{12})$\n", + "\n", + "\n", + "Avec : \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", + "$ [ \\phi_i \\phi_j | \\phi_i \\phi_j] = \\sum_a c_{ai} \\langle \\chi_a \\phi_j | \\phi_i \\phi_j \\rangle$\n", + "\n", + "$=\\sum_a \\sum_b c_{ai} c_{bj} \\langle \\chi_a \\chi_b | \\phi_i \\phi_j \\rangle $\n", + "\n", + "$=\\left(\\sum_a c_{ai} \\left(\\sum_b c_{bj} \\left(\\sum_e c_{ei} \\left(\\sum_f c_{fj} \\langle \\chi_a \\chi_b | \\chi_e \\chi_f \\rangle \\right)\\right)\\right)\\right) $\n", + "\n", + "$=\\sum_a \\sum_b \\sum_e \\sum_f \\left( c_{ai} c_{bj} c_{ei} c_{fj} \\langle \\chi_a \\chi_b | \\chi_e \\chi_f \\rangle \\right) $\n", + "\n", + "Et :\n", + "\n", + "$\\phi_i ^2 = \\left( \\sum_a c_{ai} \\chi_a \\right)^2$\n", + "$= \\sum_a c_{ai} \\sum_b c_{bi} \\chi_a \\chi_b$\n", + "\n", + "$\\phi_i ^2 - \\phi_j ^2 = \\left( \\sum_a c_{ai} \\chi_a \\right)^2 - \\left( \\sum_a c_{aj} \\chi_a \\right)^2$\n", + "$= \\sum_a c_{ai} \\sum_b c_{bi} \\chi_a \\chi_b - \\sum_a c_{aj} \\sum_b c_{bj} \\chi_a \\chi_b$\n", + "\n", + "$[\\phi_i^2 -\\phi_j^2 |\\phi_i^2 -\\phi_j^2] = [\\left( \\sum_a c_{ai} \\chi_i \\right)^2 - \\left( \\sum_a c_{aj} \\chi_j \\right)^2|\\phi_i^2 -\\phi_j^2]$\n", + "\n", + "$= \\sum_a c_{ai} \\sum_b c_{bi} [\\chi_a \\chi_b| \\phi_i^2 -\\phi_j^2 ] - \\sum_a c_{aj} \\sum_b c_{bj} [ \\chi_a \\chi_b| \\phi_i^2 -\\phi_j^2 ] $\n", + "\n", + "$= \\left(\\sum_a c_{ai} \\sum_b c_{bi} - \\sum_a c_{aj} \\sum_b c_{bj} \\right) [ \\chi_a \\chi_b| \\phi_i^2 -\\phi_j^2 ] $\n", + "\n", + "$= \\left(\\sum_a c_{ai} \\sum_b c_{bi} - \\sum_a c_{aj} \\sum_b c_{bj} \\right) [ \\chi_a \\chi_b| \\left( \\sum_a c_{ai} \\chi_i \\right)^2 - \\left( \\sum_a c_{aj} \\chi_j \\right)^2 ] $\n", + "\n", + "$= \\left(\\sum_e c_{ei} \\sum_f c_{fi} - \\sum_e c_{ej} \\sum_f c_{fj} \\right) \\left(\\sum_a c_{ai} \\sum_b c_{bi} - \\sum_a c_{aj} \\sum_b c_{bj} \\right) [ \\chi_a \\chi_b| \\chi_e \\chi_f ] $\n", + "\n", + "Mais aussi :\n", + "\n", + "$B_{ij} = [\\phi_i ^2 - \\phi_j ^2 | \\phi_i \\phi_j ] = [\\left( \\sum_a c_{ai} \\chi_i \\right)^2 - \\left( \\sum_a c_{aj} \\chi_j \\right)^2| \\phi_i \\phi_j ] $\n", + "\n", + "$= \\sum_a c_{ai} \\sum_b c_{bi} [\\chi_a \\chi_b| \\phi_i \\phi_j ] - \\sum_a c_{aj} \\sum_b c_{bj} [ \\chi_a \\chi_b| \\phi_i \\phi_j ] $\n", + "\n", + "$= \\left(\\sum_a c_{ai} \\sum_b c_{bi} - \\sum_a c_{aj} \\sum_b c_{bj} \\right) [ \\chi_a \\chi_b| \\phi_i \\phi_j ] $\n", + "\n", + "$= \\left(\\sum_a c_{ai} \\sum_b c_{bi} - \\sum_a c_{aj} \\sum_b c_{bj} \\right) \\sum_e \\sum_f c_{ei} c_{fj} [ \\chi_a \\chi_b| \\chi_e \\chi_f ] $\n", + "\n", + "Pour le calcul de $D$ (le critère de localisation qu'on maximise) , sachant que :\n", + "\n", + "$\\phi_i ^2 = \\left( \\sum_a c_{ai} \\chi_a \\right)^2$\n", + "$= \\sum_a c_{ai} \\sum_b c_{bi} \\chi_a \\chi_b$\n", + "\n", + "On aura :\n", + "\n", + "$D=\\sum_n [\\phi_n^2|\\phi_n^2]$\n", + "\n", + "$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} [\\chi_a \\chi_b| \\phi_n^2])$\n", + "\n", + "$= \\sum_n (\\sum_a c_{an} \\sum_b c_{bn} \\sum_e c_{en} \\sum_f c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$\n", + "\n", + "$= \\sum_n (\\sum_a \\sum_b \\sum_e \\sum_f c_{an} c_{bn} c_{en} c_{fn} [\\chi_a \\chi_b| \\chi_e \\chi_f])$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Localisation de Boys\n", + "\n", + "Ref:\n", + "*Localized molecular orbitals for polyatomic molecules. I. A comparison of the Edmiston‐Ruedenberg and Boys localization methods\n", + "J. Chem. Phys. 61, 3905 (1974); https://doi.org/10.1063/1.1681683\n", + "Daniel A. Kleier, Thomas A. Halgren, John H. Hall Jr., and William N. Lipscomb*\n", + "\n", + "\n", + "On procède comme pour la méthode de Edminson-Ruedenberg mais avec les intégrales $A_{12}$ et $B_{12}$ définies comme :\n", + "\n", + "$A^r_{12} = \\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n", + "$.\\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n", + "$- \\frac {1}{4}(\\langle \\phi_1 | \\bar r | \\phi_1 \\rangle $\n", + "$- \\langle \\phi_2 | \\bar r | \\phi_2 \\rangle . \\langle \\phi_1 | \\bar r | \\phi_1 \\rangle$\n", + "$- \\langle \\phi_2 | \\bar r | \\phi_2 \\rangle)$\n", + "\n", + "Et \n", + "\n", + "$B^r_{12} = (\\langle \\phi_1 | \\bar r | \\phi_1 \\rangle - \\langle \\phi_2 | \\bar r | \\phi_2 \\rangle)$\n", + "$ . \\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n", + "\n", + "Avec \n", + "\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}$\n", + "\n", + "Et le critère D à maximiser est défini tel que :\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", + "\n" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -1683,53 +1871,6 @@ "f_alpha_boys m_C;;\n" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "(* Fonction de calcul de D Boys *)\n", - "(*\n", - "let d_boys m_C = \n", - "\n", - " let phi_x_phi =\n", - " Multipole.matrix_x multipoles \n", - " |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C \n", - " in\n", - " \n", - " (*Util.debug_matrix \"phi_x_phi\" phi_x_phi;*)\n", - " \n", - " let phi_y_phi =\n", - " Multipole.matrix_y multipoles \n", - " |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C \n", - " in\n", - " \n", - " (*Util.debug_matrix \"phi_y_phi\" phi_y_phi;*)\n", - " \n", - " let phi_z_phi =\n", - " Multipole.matrix_z multipoles \n", - " |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n", - "\n", - " in\n", - " \n", - " (*Util.debug_matrix \"phi_z_phi\" phi_z_phi;*)\n", - " \n", - " let v_D_boys = \n", - " let n_mo = Mat.dim2 m_C\n", - " in\n", - " Vec.init n_mo ( fun i -> (phi_x_phi.{i,i})**2. +. (phi_y_phi.{i,i})**2. +. (phi_z_phi.{i,i})**2.)\n", - "in\n", - "Vec.sum v_D_boys;;\n", - "*)\n", - "\n", - "(*************************)\n", - "(*\n", - "Multipole.matrix_x multipoles;;\n", - "d_boys m_C;;\n", - "*)" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1874,9 +2015,11 @@ "\n", "(*************************)\n", "\n", + "(*\n", "let m_alpha,d = f_alpha m_C\n", "let alphaij = new_m_alpha m_alpha m_C 3;;\n", - "alphaij.alpha_max;;\n" + "alphaij.alpha_max;;\n", + "*)" ] }, { @@ -1913,9 +2056,11 @@ "let alpha = alphaij.alpha_max;; (* Fonction -> constante *) \n", "f_R alpha;;\n", "*)\n", + "(*\n", "alpha_v \"deloc\" alphaij;;\n", "let alpha = (alpha_v \"loc\" alphaij) ;;\n", - "f_R alpha ;;\n" + "f_R alpha ;;\n", + "*)" ] }, { @@ -2184,7 +2329,7 @@ " let alpha = (alpha_v loc_deloc alphaij) *. epsilon (* Fonction -> constante *)\n", " in\n", "\n", - " Printf.printf \"%f\\n%!\" alpha;\n", + " (*Printf.printf \"%f\\n%!\" alpha;*)\n", " \n", " (* Indice i et j du alpha max après calcul *)\n", " let indice_i = alphaij.indice_ii (* Fonction -> constante *)\n", @@ -2252,7 +2397,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "jupyter": { + "outputs_hidden": true + } + }, "outputs": [], "source": [ "(* Calcul *)\n", @@ -2264,7 +2413,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "jupyter": { + "outputs_hidden": true + } + }, "outputs": [], "source": [ "final_m_C m_C \"boys\" \"loc\" 1. 1 0. 10e-7;;" @@ -2273,7 +2426,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "jupyter": { + "outputs_hidden": true + } + }, "outputs": [], "source": [ "final_m_C new_m \"Boys\" \"deloc\" 1. 1 0. 10e-7;;\n" @@ -2285,19 +2442,17 @@ "metadata": {}, "outputs": [], "source": [ - "nocc;;\n", - "let toto = Vec.init (nocc) (fun i -> float_of_int(i));;\n", - "let vec_list = Vec.to_list toto;;\n", - "let g a = int_of_float(a);;\n", - "let tutu = List.map g vec_list\n", - "\n", - "let int_list vec = \n", - " let float_list = Vec.to_list vec\n", - " in\n", - " let g a = int_of_float(a)\n", - "in List.map g float_list;;\n", - "\n", - "int_list toto;;" + "let m_occ , m_vir = split_mat new_m list_occ;;" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "let loc_m_occ = final_m_C m_occ \"boys\" \"loc\" 1. 1 0. 10e-7;;\n", + "let loc_m_vir = final_m_C m_vir \"boys\" \"loc\" 1. 1 0. 10e-7;;" ] }, {