22 KiB
Rotation d'orbitales¶
On peut écrire une matrice de rotation $R$ comme $$ R = \exp(X) $$ où X est une matrice réelle anti-symétrique ($X_{ij} = -X_{ji}$) où $X_{ij}$ est la valeurs de l'angle de la rotation entre $i$ et $j$. Par exemple: $$ R = \left( \begin{array}{cc} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array} \right) = \exp \left( \begin{array}{cc} 0 & \theta \
- \theta & 0 \end{array} \right) $$
D'abord on diagonalise $X^2$ : $$ X^2 = W(-\tau^2) W^\dagger. $$ Ensuite, $R$ peut être écrit comme $$ R = W \cos(\tau) W^\dagger + W \tau^{-1} \sin (\tau) W^\dagger X $$
TODO¶
Ecrire une fonction qui calcule $R$ a partir de $X$. Tu vas avoir besoin de Lacaml:
- Intro a Lapack : https://lipn.univ-paris13.fr/~coti/doc/tutolapack.pdf
- Lapack user's guide : http://www.netlib.org/lapack/lug/
- Lacaml : https://mmottl.github.io/lacaml/
Pour diagonaliser une matrice : syev
Pour multiplier des matrices : gemm
Localisation de Boys¶
$N$ orthonormal molecular orbitals $$ \phi_i({\bf r})= \sum_{k=1}^N c_{ik} \chi_k $$ $$ {\cal B}_{2|4} = \sum_{i=1}^N \langle \phi_i | ( {\bf \vec{r}} - \langle \phi_i| {\bf \vec{r}} | \phi_i …\rangle)^{2|4} | \phi_i \rangle $$ $$ {\cal B}_2= \sum_{i=1}^N \big[ \langle x^2 \rangle_i - \langle x \rangle^2_i + \langle y^2 \rangle_i - …\langle y \rangle^2_i + \langle z^2 \rangle_i - \langle z \rangle^2_i \big] $$ $$ {\cal B}4= \sum{i=1}^N \big[ \langle x^4 \rangle_i - 4 \langle x^3 \rangle_i \langle x \rangle_i
- 6 \langle x^2 \rangle_i \langle x \rangle^2_i
- 3 \langle x \rangle^4_i \big] + \big[ ...y...] + \big[ ...z...] $$ Minimization of ${\cal B}$ with respect to an arbitrary rotation $R$ $$ \langle R \phi_i x^n R \phii \rangle = \sum{k,l=1}^N R{ik} R{il} \langle \phi_k| x^n | \phil …\rangle= $$ $$ \sum{k,l=1}^N R{ik} R{il} \sum{m,o=1}^N c{km} c_{ln} \langle \chi_m | x^n |\chio \rangle $$ We need to compute $$ S^x{n;mo}= \langle \chi_m | x^n |\chi_o \rangle $$
Rotation de deux orbitales¶
Rotation of angle $\theta$ $$ \tilde{\phi}_1 = cos\theta \phi_1 -sin\theta \phi_2 $$ $$ \tilde{\phi}_2 = sin\theta \phi_1 +cos\theta \phi_2 $$ Let us note $$ {\cal B}_x(\theta) = \langle x^2 \rangle_{\tilde{1}} - \langle x \rangle^2_{\tilde{1}} + \langle x^2 \rangle_{\tilde{2}} - \langle x \rangle^2_{\tilde{2}} $$ and \begin{eqnarray} A^x_{ij} & = & \langle \phi_i| x^2 | \phi_j \rangle \\ B^x_{ij} & = & \langle \phi_i| x | \phi_j \rangle \;\; i,j=1,2 \end{eqnarray} We have $$ {\cal B}_x(\theta) = A^x_{11}+A^x_{22} - ({\tilde B}^{x2}_{11} + {\tilde B}^{x2}_{22}) $$ with $$ {\tilde B}^{x2}_{11}= (cos^2\theta B^x_{11} + sin^2\theta B^x_{22} - sin2\theta B^x_{12})^2 $$ $$ {\tilde B}^{x2}_{22}= (sin^2\theta B^x_{11} + cos^2\theta B^x_{22} + sin2\theta B^x_{12})^2 $$ $$ {\cal B}(\theta)= {\cal B}_x(\theta)+{\cal B}_y(\theta)+{\cal B}_z(\theta)=(x) + (y)+(z) $$ with $$ {\cal B}x(\theta)= A^x{11}+A^x_{22}
- [(cos^4\theta+ sin^4\theta)({B^x{11}}^2+ {B^x{22}}^2 ) $$ $$
- 2 sin^2 2\theta {B^x_{12}}^2
- 2 sin 2\theta cos 2\theta (({B^x{22}} -{B^x{11}} ) {B^x_{12}}] $$ and idem for $y$ and $z$. Using the fact that $$ cos^4\theta+ sin^4\theta= \frac{1}{4} ( 3 + cos4\theta) $$ $$ {\cal B}x(\theta)= A^x{11}+A^x_{22}
- [ \frac{1}{4} ( 3 + cos4\theta)({B^x{11}}^2+ {B^x{22}}^2 ) $$ $$
- (1 -cos 4\theta) {B^x_{12}}^2
- sin 4\theta (({B^x{22}} -{B^x{11}} ) {B^x{12}}] $$ Finally, we get \begin{equation} {\cal B}(\theta)= {\cal B}(0) + \frac{1}{4} [(1-cos4\theta)\beta+sin4\theta \gamma] \end{equation} where $$ {\cal B}(0)= A^x{11}+A^x{22} -((B^{x}{11})^2+(B^{x}{22})^2) + [...y...] + [...z...] $$ $$ \beta= (B^x{11}-B^x{22})^2 - 4 {(B^x{12})}^2 + [...y...] + [...z...] $$ and $$ \gamma= 4 B^x{12} (B^{x}{11}-B^x_{22}) + [...y...] + [...z...] $$ Let us compute the derivative; we get $$ \frac{\partial {\cal B}(\theta)}{\partial \theta} = \beta sin4\theta
- \gamma cos4\theta $$ Extrema of ${\cal B}(\theta)$ \begin{equation} tg4\theta= -\frac{\gamma}{\beta} \end{equation} There are four extrema: $$ 4\theta; \;\; 4\theta +\pi; \;\; 4\theta+ 2\pi; \;\; 4\theta+ 3\pi $$ Value of the second derivative of $\cal{B}$ at the extrema Value of $\cal{B}$ at the extrema \begin{equation} \frac{\partial^2 B(\theta)}{\partial \theta^2}= 4 cos4\theta \frac{\beta^2 + \gamma^2}{\beta} \end{equation} There are two minima and two maxima since $cos4\theta= -cos4(\theta+\pi)= -cos4(\theta+2\pi)=-cos4(\theta+3\pi)$. Value of $\cal{B}$ at the extrema \begin{equation} {\cal B}(\theta)= {\cal B}(0) + \frac{1}{4} (\beta -\frac{\beta^2 + \gamma^2}{\beta} {cos4\theta}) \end{equation}
#use "topfind";;
#require "jupyter.notebook";;
#require "lacaml.top";;
#require "alcotest";;
#require "str";;
#require "bigarray";;
#require "zarith";;
#require "getopt";;
open Lacaml.D
Construction de la matrice de rotation R de taille n par n:
On part de X ( n * n ) une matrice antisymétrique avec les éléments diagonaux nuls et où le terme X_ij correspond à la valeur de l'angle de rotation i et j (dans notre cas entre l'orbitale i et j).
Tout d'abord on met au carré cette matrice X ( gemm X X ) tel que :
X² = X . X
Ensuite on la diagonalise la matrice X² ( syev X²). Après la diagonalidsation on obtient une matrice contenant les vecteurs propres ( W ) et un vecteur contenant les différentes valeurs propres. Ces valeurs propres correspondent aux valeurs propres de X² et sont négatives, pour obtenir celle correpondant à X et qui sont positives on doit appliquer la fonction valeur absolue au vecteur et ensuite prendre sa racine carré. On obtient donc un nouveau vecteur contenant les valeurs propres de X. Il faut ensuite transformer ce vecteur en une matrice diagonale, pour cela on construit une matrice (tau) n par n où les éléments extra diagonaux sont nulles et les éléments diagonaux prennent les différentes valeurs du vecteur, tel que l'élément 1 du vecteur devient l'élément {1,1} de la matrice.
On peut ensuite écrire R comme :
R = W cos(tau) W_t + W tau⁻¹ sin(tau) W_t X(*
Construction de la matrice de rotation R pour la rotation de n orbitales
*)
(*
NB: il faut que je change les noms des matrices car ce n'est pas toujours clair
*)
(* Construction de la matrice X n par n *)
let n =2
let ran=Mat.random ~range:1. ~from:0. n n;; (* Construction d'une matrice random n*n *)
(* Antisymétrisation de la matrice *)
let antisym = Mat.init_cols n n (fun i j -> if i>j then -. ran.{i,j} else ran.{j,i});;
(* 0 sur la diagonale -> X *)
let m = Mat.init_cols n n (fun i j -> if i=j then 0. else antisym.{i,j});;
(* Fonction 2 en un *)
(*
let x = Mat.init_cols n n (fun i j -> if i=j then 0.
else if i>j then -. ran.{i,j} else ran.{j,i});;
*)
(* Mise au carré de la matrice : X -> X² *)
let mm = gemm m m;;
(* Copie de mm qui va être écrasée : X² *)
let copie_mm = lacpy mm;;
(* Diagonalisation de X² *)
let diag = syev mm;;
copie_mm;; (* Copie de la matrice avant diagonalisation *)
mm;; (* Matrice avec les eigenvectors -> W *)
diag;; (* Vecteur contenant les eigenvalues -> -tau² *)
(* Passage de -tau² à tau² *)
let square_tau= Vec.abs diag;;
(* Passage de tau² à tau *)
let tau = Vec.sqrt square_tau;;
(* Passage du vecteur tau à la matrice avec les valeurs propres sur la diagonale *)
let m_tau = Mat.init_cols n n (fun i j -> if i=j then tau.{i} else 0.);;
(*Test opération sur tau, fonction générale*)
let op_tau f = Mat.init_cols n n (fun i j -> if i=j then f tau.{i} else 0.);;
let m_cos_tau = op_tau cos ;;
let m_sin_tau = op_tau sin ;;
(* Calcul de cos tau à partir du vecteur tau *)
let cos_tau = Vec.cos tau;;
(* Matrice cos tau *)
let m_cos_tau = Mat.init_cols n n (fun i j -> if i=j then cos_tau.{i} else 0.)
(* Calcul de sin tau à partir du vecteur tau *)
let sin_tau = Vec.sin tau;;
(* Matrice de sin tau *)
let m_sin_tau = Mat.init_cols n n (fun i j -> if i=j then sin_tau.{i} else 0.)
(* Calcul de la transposée de W (mm) *)
let transpose_mm = Mat.transpose_copy mm;;
mm;;(* Vérification *) (* OK *)
(* Calcul du vecteur tau^-1 *)
let tau_1 = Vec.init n (fun i -> 1. /. tau.{i}) ;;
tau;; (* Vétification *) (* OK *)
(* Calcul de la matrice tau^⁻1 *)
let m_tau_1 = Mat.init_cols n n (fun i j -> if i=j then tau_1.{i} else 0.);;
(* Calcul de R *)
(*
let r = gemm mm (gemm m_cos_tau transpose_mm) + gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m)))
*)
(* gemm mm (gemm m_cos_tau transpose_mm) -> W cos(tau) Wt *)
let a = gemm mm (gemm m_cos_tau transpose_mm);;
(* gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m))) -> W tau^-1 sin(tau) Wt X *)
let b = gemm mm (gemm m_tau_1 ( gemm m_sin_tau (gemm transpose_mm m)));;
(* Somme de 2 matrices, ici -> a + b -> R *)
let m_r = Mat.add a b;; (* Matrice de rotation R *)
gesvd m_r;;
Localisation des orbitales par rotation de paires d'orbitales
Ref : Edmiston, C., & Ruedenberg, K. (1963). Localized Atomic and Molecular Orbitals. Reviews of Modern Physics, 35(3), 457–464. doi:10.1103/revmodphys.35.457
Méthode :
Le but de cette méthode est de maximiser D :
D(phi) = Somme sur n de < Phi²_n | 1/r_12 | Phi²_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.
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 :
i~ (x) = cos(gamma) i(x) + sin(gamma) j(x)
j~ (x) = - sin(gamma) i(x) - cos(gamma) j(x)
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 :
cos(4 alpha) = -A_12/(A²_12 + B²_12)^(1/2)
<=> alpha = (1/4) acos {-A_12/(A²_12 + B²_12)^(1/2)}
sin(4 alpha) = B_12/(A²_12 + B²_12)^(1/2)
<=> alpha = (1/4) asin {B_12/(A²_12 + B²_12)^(1/2)}
tan(4 alpha) = -B_12/A_12
<=> alpha = (1/4) atan {-B_12/A_12}
Avec : $A_{ij} = [ \phi_i \phi_j | \phi_i \phi_j] - 1/4 [\phi_i^2 - \phi_j^2 | \phi_i^2 - \phi_j^2] $
Où :
$ [ \phi_i \phi_j | \phi_i \phi_j] = \sum_a c_{ai} \langle \chi_a \phi_j | \phi_i \phi_j \rangle$
$=\sum_a \sum_b c_{ai} c_{bj} \langle \chi_a \chi_b | \phi_i \phi_j \rangle $
$=\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) $
$=\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) $
Et :
$\phi_i ^2 = \left( \sum_a c_{ai} \chi_a \right)^2$ $= \sum_a c_{ai} \sum_b c_{bi} \chi_a \chi_b$
$\phi_i ^2 - \phi_j ^2 = \left( \sum_a c_{ai} \chi_a \right)^2 - \left( \sum_a c_{aj} \chi_a \right)^2$ $= \sum_a c_{ai} \sum_b c_{bi} \chi_a \chi_b - \sum_a c_{aj} \sum_b c_{bj} \chi_a \chi_b$
$[\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]$
$= \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 ] $
$= \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 ] $
$= \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 ] $
$= \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 ] $
Mais aussi :
$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 ] $
$= \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 ] $
$= \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 ] $
$= \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 ] $
On extrait de la matrice des OMs les orbitales i et j (vecteurs colonnes) pour former une nouvelle matrice Ksi de dimension n * 2.
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 :
R = ( cos(alpha) -sin(alpha) ) ( sin(alpha) cos(alpha) )
On applique R à Ksi : R Ksi = Ksi~
On obtient Ksi~ la matrice contenant les deux nouvelles OMs i~ et j~ obtenues par rotation de i et j.
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~. Pour cela on créer des matrices intérmédiaires:
- une matrice (m_ksi) n par n avec tous ses éléments nuls sauf les colonnes qui contiennent les OMs i et j
- une matrice (m_ksi_tilde) n par n avec tous ses éléments nuls sauf les colonnes qui contiennent les OMs i~ et j~ Ainsi, on soustrait à la matrice des OMs Phi la matrice m_ksi pour supprimer les OMs i et j de celle ci puis on additionne cette nouvelle matrice à la matrice m_ksi_tilde pour créer la nouvelle matrice des OMs Phi~.
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, c'est à dire maximiser le terme de répulsion.
let a= 1.;;
(* Matrice de rotation 2*2 *)
(*let rot = Mat.init_cols 2 2 (fun i j -> match (i,j) with
| (1,1) -> cos a
| (2,2) -> cos a
| (2,1) -> sin a
| (1,2) -> -. sin a );;*)
let rot = Mat.init_cols 2 2 (fun i j -> if i=j then cos a
else if i>j then sin a
else -. sin a );;
(* Fonction d'extraction des 2 vecteurs propres dans la matrice (Phi~) n par 2 pour appliquer R afin d'effectuer
la rotation des orbitales *)
ran.{1,2};; (*{1,2} -> 1ere ligne, 2e colonne*)
let p = 1;;
let q = 2;;
let ksi = Mat.init_cols n 2 (fun i j -> if j=1 then ran.{i,p} else ran.{i,q} );;
(* Fonction de calcul de ksi~ *)
let ksi_tilde = gemm ksi rot;;
(* Pour la réinjection on créer des matrices intérmédiares, une matrice nulle partout sauf sur
les colonnes de i et j et de i~ et j~. On fait la différence de la première matrice avec la matrice
des OMs Phi afin de substituer les colonnes de i et j par des zéro et ensuite sommer cette matrice avec
celle contenant i~ et j~ *)
(* Matrice intérmédiare pour l'injection de ksi~ *)
let m_ksi_tilde = Mat.init_cols n n (fun i j -> if j=q then ksi_tilde.{i,p}
else if j=p then ksi_tilde.{i,q}
else 0.);;
(* Matrice intermédiaire pour supprimer ksi *)
let m_ksi = Mat.init_cols n n (fun i j -> if j=q then ksi.{i,p}
else if j=p then ksi.{i,q}
else 0.);;
(* Test soustraction *)
let test = Mat.sub ran m_ksi;;
(* Construction de la nouvelle matrice d'OMs phi~ *)
let new_ran = Mat.add m_ksi_tilde test;;
let toto = Mat.max2 new_ran
let vec_test = Mat.as_vec new_ran
|> iamax;;
11 mod 5;;