scemama/StageYann
1
0
Fork 0
Code Pull Requests Releases Activity

Boyls : rotation de paires

Browse Source
This commit is contained in:
Yann Damour 2020-04-24 17:07:06 +02:00
parent 55a1a68297
commit 5c48dc4c6c
1 changed files with 149 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

151
Boys.ipynb
View File

@ -202,6 +202,8 @@
"metadata": {},
"outputs": [],
"source": [
"#use \"topfind\";;\n",
"#require \"jupyter.notebook\";;\n",
"#require \"lacaml.top\";;\n",
"#require \"alcotest\";;\n",
"#require \"str\";;\n",
@ -217,6 +219,9 @@
"metadata": {},
"outputs": [],
"source": [
"(*\n",
"Construction de la matrice de rotation R pour la rotation de n orbitales \n",
"*)\n",
"(*\n",
"NB: il faut que je change les noms des matrices car ce n'est pas toujours clair\n",
"*)\n",
@ -328,8 +333,150 @@
"(* gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) -> W tau^-1 sin(tau) Wt X *)\n",
"let b = gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m)));;\n",
"\n",
"(* Fonction somme de 2 matrices, ici -> a + b -> R *)\n",
"let sum_m = Mat.init_cols n n (fun i j -> a.{i,j} +. b.{i,j});;\n"
"(* Somme de 2 matrices, ici -> a + b -> R *)\n",
"let m_r = Mat.add a b;; (* Matrice de rotation R *)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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",
"\n",
"Méthode :\n",
"\n",
"On va créer des nouvelles orbitales i~ et 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",
" j~ (x) = sin(gamma) i(x) - cos(gamma) j(x)\n",
"\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 :\n",
"\n",
" cos(4 alpha) = -A_12/(A²_12 + B²_12)^(1/2)\n",
" <=> alpha = (1/4) acos {-A_12/(A²_12 + B²_12)^(1/2)}\n",
"\n",
" sin(4 alpha) = B_12/(A²_12 + B²_12)^(1/2)\n",
" <=> alpha = (1/4) asin {B_12/(A²_12 + B²_12)^(1/2)}\n",
"\n",
" tan(4 alpha) = -B_12/A_12\n",
" <=> alpha = (1/4) atan {-B_12/A_12}\n",
"\n",
"Avec : \n",
" A_12 = ...\n",
" B_12 = ...\n",
"/ *\n",
"Dans la démonstration de Michel les définitions de B(theta), A_ij et B_ij sont claires, mais j'ai du mal avec les définitions utilisées par Edmiston, C., & Ruedenberg, K. Je pense que je comprend mal la notation car je ne vois pas comment les retrouver en partant de la définition de D(phi) et de la définition des nouvelles orbitales créées par la rotation des anciennes.\n",
"\n",
"D(phi) est défini tel que : \n",
" \n",
" D(phi) = Somme sur n de [ phi²_n | phi²_n ]\n",
"\n",
"Le but est de trouvé une transformation :\n",
" \n",
" λ_v(x) = Somme sur n de phi_n(x) T_nv\n",
" \n",
"tel que : \n",
"\n",
" D(λ) = Somme sur v de [λ²_v|λ²_v]\n",
" \n",
"Afin de maximiser D.\n",
"\n",
"Cas 2D:\n",
"\n",
" i~ (x) = cos(gamma) i(x) - sin(gamma) j(x)\n",
" j~ (x) = sin(gamma) i(x) - cos(gamma) j(x)\n",
"\n",
" ? => ? D(~) = D(gamma) + A_12 + ( A²_12 + B²_12 )^(1/2) cos 4(gamma-alpha)\n",
"\n",
"\n",
"* /\n",
"\n",
"\n",
"\n",
"\n",
"On extrait de la matrice des OMs les orbitales i et j (vecteurs colonnes) pour former une nouvelle matrice Ksi de dimension n * 2. \n",
"\n",
"Pour effectuer la rotation des orbitales i et j on utilise la matrice de rotation R pour la rotation de 2 orbitales, qui est définie comme :\n",
"\n",
"R = ( cos(alpha) -sin(aplha) )\n",
" ( sin(alpha) cos(alpha) )\n",
" \n",
"On applique R à Ksi : R Ksi = Ksi~ \n",
"\n",
"On obtient Ksi~ la matrice contenant les deux nouvelles OMs i~ et j~ obtenues par rotation de i et j.\n",
"\n",
"On réinjecte ces deux nouvelles orbitales i~ et j~ à la place des anciennes orbitales i et j dans la matrice des OMs Phi, ce qui nous donne une nouvelle matrice des OMs Phi~.\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. (Critère de convergence ? sur alpha ?)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let a= 1.;;\n",
"\n",
"(* Matrice de rotation 2*2 *)\n",
"(*let rot = Mat.init_cols 2 2 (fun i j -> match (i,j) with\n",
" | (1,1) -> cos a\n",
" | (2,2) -> cos a\n",
" | (2,1) -> sin a\n",
" | (1,2) -> -. sin a );;*)\n",
" \n",
"let rot = Mat.init_cols 2 2 (fun i j -> if i=j then cos a\n",
" else if i>j then sin a \n",
" else -. sin a );;\n",
"\n",
"(* Fonction *)\n",
"ran.{1,2};; (*{1,2} -> 1ere ligne, 2e colonne*)\n",
"let p = 1;;\n",
"let q = 2;;\n",
"let phi_tilde = Mat.init_cols n 2 (fun i j -> if j=1 then ran.{i,p} else ran.{i,q} );;"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Fonction de calcul de ksi~ *)\n",
"let ksi_tilde = gemm ksi rot;;"
]
},
{
"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~ *)\n",
"let m_ksi_tilde = Mat.init_cols n n (fun i j -> if j=q then ksi_tilde.{i,p}\n",
" else if j=p then ksi_tilde.{i,q}\n",
" else 0.);;\n",
"(* Matrice intermédiaire pour supprimer ksi *) \n",
"let m_ksi = Mat.init_cols n n (fun i j -> if j=q then ksi.{i,p}\n",
" else if j=p then ksi.{i,q}\n",
" else 0.);;\n",
"(* Test soustraction -> problème de type sur new_ran après*)\n",
"let test = Mat.sub ran m_ksi;;\n",
"\n",
"(*Construction de la nouvelle matrice d'OMs phi~*)\n",
"let new_ran = Mat.sum m_ksi_tilde test;; "
]
}
],