46 KiB
Installation de QCaml¶
- Clonage de QCaml:
git clone https://gitlab.com/scemama/QCaml.git
- Installation des dependances:
opam install ocamlbuild ocamlfind lacaml gnuplot getopt alcotest zarith cd QCaml ./configure
- Compilation
make
Liens utiles¶
- Documentation des modules : https://sanette.github.io/ocaml-api/4.10/index.html
- Lacaml : http://mmottl.github.io/lacaml/api/lacaml/
- Gnuplot : https://github.com/c-cube/ocaml-gnuplot
Initialisation¶
Bloc a executer avant de pouvoir utiliser QCaml dans le Notebook
let png_image = print_endline ;;
(* --------- *)
(*Mettre le bon chemin ici *)
#cd "/home/ydamour/QCaml";;
#use "topfind";;
#require "jupyter.notebook";;
#require "gnuplot";;
let png_image name =
Jupyter_notebook.display_file ~base64:true "image/png" ("Notebooks/images/"^name)
;;
#require "lacaml.top";;
#require "alcotest";;
#require "str";;
#require "bigarray";;
#require "zarith";;
#require "getopt";;
#directory "_build";;
#directory "_build/Basis";;
#directory "_build/CI";;
#directory "_build/MOBasis";;
#directory "_build/Nuclei";;
#directory "_build/Parallel";;
#directory "_build/Perturbation";;
#directory "_build/SCF";;
#directory "_build/Utils";;
#load "Constants.cmo";;
#load_rec "Util.cma";;
#load_rec "Matrix.cmo";;
#load_rec "Simulation.cmo";;
#load_rec "Spindeterminant.cmo";;
#load_rec "Determinant.cmo";;
#load_rec "HartreeFock.cmo";;
#load_rec "MOBasis.cmo";;
#load_rec "F12CI.cmo";;
#install_printer AngularMomentum.pp_string ;;
#install_printer Basis.pp ;;
#install_printer Charge.pp ;;
#install_printer Coordinate.pp ;;
#install_printer Vector.pp;;
#install_printer Matrix.pp;;
#install_printer Util.pp_float_2darray;;
#install_printer Util.pp_float_array;;
#install_printer Util.pp_matrix;;
#install_printer HartreeFock.pp ;;
#install_printer Fock.pp ;;
#install_printer MOClass.pp ;;
#install_printer DeterminantSpace.pp;;
#install_printer SpindeterminantSpace.pp;;
#install_printer Bitstring.pp;;
let pp_mo ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;
#install_printer pp_mo;;
(* --------- *)
open Lacaml.D
Calculs¶
H$_2$O¶
let xyz_string =
"3
Water
O 0. 0. 0.
H -0.756950272703377558 0. -0.585882234512562827
H 0.756950272703377558 0. -0.585882234512562827
"
Base atomique¶
Les bases atomiques sont telechargeables ici: https://www.basissetexchange.org/
On telecharge la base cc-pvdz depuis le site:
let basis_string = "
HYDROGEN
S 4
1 1.301000E+01 1.968500E-02
2 1.962000E+00 1.379770E-01
3 4.446000E-01 4.781480E-01
4 1.220000E-01 5.012400E-01
S 1
1 1.220000E-01 1.000000E+00
P 1
1 7.270000E-01 1.0000000
OXYGEN
S 9
1 1.172000E+04 7.100000E-04
2 1.759000E+03 5.470000E-03
3 4.008000E+02 2.783700E-02
4 1.137000E+02 1.048000E-01
5 3.703000E+01 2.830620E-01
6 1.327000E+01 4.487190E-01
7 5.025000E+00 2.709520E-01
8 1.013000E+00 1.545800E-02
9 3.023000E-01 -2.585000E-03
S 9
1 1.172000E+04 -1.600000E-04
2 1.759000E+03 -1.263000E-03
3 4.008000E+02 -6.267000E-03
4 1.137000E+02 -2.571600E-02
5 3.703000E+01 -7.092400E-02
6 1.327000E+01 -1.654110E-01
7 5.025000E+00 -1.169550E-01
8 1.013000E+00 5.573680E-01
9 3.023000E-01 5.727590E-01
S 1
1 3.023000E-01 1.000000E+00
P 4
1 1.770000E+01 4.301800E-02
2 3.854000E+00 2.289130E-01
3 1.046000E+00 5.087280E-01
4 2.753000E-01 4.605310E-01
P 1
1 2.753000E-01 1.000000E+00
D 1
1 1.185000E+00 1.0000000
"
Une orbitale atomique centree sur l'atome A est composee d'une contraction de m Gaussiennes: $$ \chi(r) = \sum_{i=1}^m a_i \exp \left( -\alpha_i |r-r_A|^2 \right) $$ Dans le fichier de base, la 2ieme colonne represente l'exposant $\alpha_i$ et la 3ieme colonne le coefficient de contraction $a_i$.
let nuclei =
Nuclei.of_xyz_string xyz_string
let basis =
Basis.of_nuclei_and_basis_string nuclei basis_string
let simulation =
Simulation.make ~charge:0 ~multiplicity:1 ~nuclei basis
let ao_basis =
Simulation.ao_basis simulation
Plot des orbitales atomiques¶
let plot data filename =
let output = Gnuplot.Output.create (`Png filename) in
let gp = Gnuplot.create () in
Gnuplot.set gp ~use_grid:true;
List.map Gnuplot.Series.lines_xy data
|> Gnuplot.plot_many gp ~output;
Gnuplot.close gp ;
Jupyter_notebook.display_file ~base64:true "image/png" filename
;;
let x_values =
let n = 1000 in
let xmin, xmax =
let coord =
Array.map snd nuclei
|> Array.map (fun a -> Coordinate.(get X) a)
in
Array.sort compare coord;
coord.(0) -. 4. ,
coord.(Array.length coord -1) +. 4.
in
let dx =
(xmax -. xmin) /. (float_of_int n -. 1.)
in
Array.init n (fun i -> xmin +. (float_of_int i)*.dx)
in
let data =
Array.map (fun x ->
let point = Coordinate.make_angstrom
{ Coordinate.
x ; y = 0. ; z = 0.
} in
AOBasis.values ao_basis point
) x_values
|> Mat.of_col_vecs
|> Mat.transpose_copy
|> Mat.to_col_vecs
|> Array.map Vec.to_list
|> Array.map (fun l -> List.mapi (fun i y -> (x_values.(i),y)) l)
|> Array.to_list
in
plot data "test_data.png"
Calcul Hartree-Fock¶
let hf = HartreeFock.make simulation
let mo_basis = MOBasis.of_hartree_fock hf
Orbitales moleculaires : $$ \phi_j(r) = \sum_{i=1}^{N_b} C_{ij} \chi_i(r) $$
- $i$: lignes
- $j$: colonnes
Extraction des OM de la structure de donnees mo_basis comme une matrice $C$ utilisable avec Lacaml:
let mo_coef = MOBasis.mo_coef mo_basis
let plot data filename =
let output = Gnuplot.Output.create (`Png filename) in
let gp = Gnuplot.create () in
Gnuplot.set gp ~use_grid:true;
List.map Gnuplot.Series.lines_xy data
|> Gnuplot.plot_many gp ~output;
Gnuplot.close gp ;
Jupyter_notebook.display_file ~base64:true "image/png" filename
;;
let x_values =
let n = 1000 in
let xmin, xmax =
let coord =
Array.map snd nuclei
|> Array.map (fun a -> Coordinate.(get X) a)
in
Array.sort compare coord;
coord.(0) -. 4. ,
coord.(Array.length coord -1) +. 4.
in
let dx =
(xmax -. xmin) /. (float_of_int n -. 1.)
in
Array.init n (fun i -> xmin +. (float_of_int i)*.dx)
in
let data =
let result =
Array.map (fun x ->
let point = Coordinate.make_angstrom
{ Coordinate.
x ; y = 0. ; z = 0.
} in
MOBasis.values mo_basis point
) x_values
|> Mat.of_col_vecs
|> Mat.transpose_copy
|> Mat.to_col_vecs
|> Array.map Vec.to_list
|> Array.map (fun l -> List.mapi (fun i y -> (x_values.(i),y)) l)
in
[ result.(0) ; result.(1) ; result.(2) ]
in
plot data "test_data.png"
Calcul Hartree-Fock a la main¶
Methode¶
$$ \hat{H} = \hat{T} + \hat{V}^{\text{NN}} + \hat{V}^{\text{eN}} + \hat{V}^{\text{ee}} $$On exprime les differentes quantites dans la base des orbitales atomiques pour chacun des termes:
\begin{eqnarray} T_{ij} & = & \langle \chi_i | \hat{T} | \chi_j \rangle = \iiint \chi_i(\mathbf{r}) \left( -\frac{1}{2} \Delta \chi_j(\mathbf{r}) \right) d\mathbf{r} \\ V^{NN} & = & \text{constante} \\ V^{eN}_{ij} & = & \langle \chi_i | \hat{V}^{\text{eN}} | \chi_j \rangle = \sum_A \iiint \chi_i(\mathbf{r}) \frac{-Z_A}{|\mathbf{r} - \mathbf{R}_A|} \chi_j(\mathbf{r}) d\mathbf{r} \\ V^{ee}_{ijkl} & = & \langle \chi_i \chi_j | \hat{V}^{\text{ee}} | \chi_k \chi_l \rangle = \iiint \iiint \chi_i(\mathbf{r}_1) \chi_j(\mathbf{r}_2) \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} \chi_k(\mathbf{r}) \chi_l(\mathbf{r}) d\mathbf{r}_1 d\mathbf{r}_2 \\ \end{eqnarray}ou $\mathbf{R}_A$ est la position du noyau $A$ et $Z_A$ sa charge.
Dans la methode Hartree-Fock, l'electron 1 ne "voit" pas directement l'electron 2, mais il voit l'electron 2 comme une densite de charge.
On definit donc 2 operateurs pour chaque orbitale:
- Coulomb : $$ \hat{J}_j \chi_i(\mathbf{r}_1) = \chi_i(\mathbf{r}_1) \iiint \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} |\chi_j(\mathbf{r}_2)|^2 d \mathbf{r}_2 $$
- Echange : $$ \hat{K}_j \chi_i(\mathbf{r}_1) = \chi_j(\mathbf{r}_1) \iiint \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} \chi_i(\mathbf{r}_2) \chi_j(\mathbf{r}_2) d \mathbf{r}_2 $$ et on n'a plus que des operateurs a 1 electron.
ou $P$ est la matrice densite definie comme $$ P_{ij} = \sum_{k=1}^{N_{\text{occ}}} 2 C_{ik} C_{kj} $$ et $C$ est la matrice des coefficients des orbitales moleculaires exprimees dans la base des orbitales atomiques.
Une orbitale moleculaire est une combinaison lineaire d'orbitale atomique: $$ \phi_k(\mathbf{r}) = \sum_i C_{ik} \chi_i(\mathbf{r}) $$
Les orbitales moleculaires sont toutes orthogonales entre elles, et normalisees a 1.
La methode Hartree-Fock permet d'obtenir les orbitales moleculaires qui minimisent l'energie quand la fonction d'onde est exprimee comme un determinant de Slater: $$ \Psi(\mathbf{r}_1,\dots,\mathbf{r}_N) = \left| \begin{array}{ccc} \phi_1(\mathbf{r}_1) & \dots & \phi_N(\mathbf{r}_1) \\ \vdots & \ddots & \vdots \\ \phi_1(\mathbf{r}_N) & \dots & \phi_N(\mathbf{r}_N) \end{array}\right| $$
Application¶
On commence par creer une base orthogonale $X$ a partir des OA
let m_X =
AOBasis.ortho ao_basis
Pour pouvoir calculer les matrices $J$ et $K$, il faut avoir une matrice densite. On part donc avec une matrice densite d'essai. Une facon de faire est de diagonaliser l'Hamiltonien sans $V_{ee}$.
let nocc = 5 (* On a 10 electrons, donc 5 orbitales occupees. *)
let c_of_h m_H =
(* On exprime H dans la base orthonormale *)
let m_Hmo =
Util.xt_o_x m_H m_X (* H_mo = X^t H X *)
in
(* On diagonalise cet Hamiltonien *)
let m_C', _ =
Util.diagonalize_symm m_Hmo
in
(* On re-exprime les MOs dans la base des AOs (non-orthonormales) *)
gemm m_X m_C' (* C = X.C' *)
let m_C =
match Guess.make ~nocc ~guess:`Hcore ao_basis with
| Hcore m_H -> c_of_h m_H
| _ -> assert false
On construit la matrice densite
let m_P =
(* P = 2 C.C^t *)
gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C
On construit les matrices $H_c = T+V_{\text{eN}}$, et $J$ et $K$:
let m_Hc, m_J, m_K =
let f =
Fock.make_rhf ~density:m_P ao_basis
in
Fock.(core f, coulomb f, exchange f)
On construit l'operateur de Fock: $ F = H_c + J - K $
let m_F =
Mat.add m_Hc (Mat.sub m_J m_K)
On l'exprime dans la base orthonormale, on le diagonalise, et on le re-exprime ses vecteurs propres dans la base d'OA.
let m_C =
let m_C', _ =
Util.xt_o_x m_F m_X
|> Util.diagonalize_symm
in
gemm m_X m_C'
On calcule l'energie avec ces nouvelles orbitales:
let energy =
(Simulation.nuclear_repulsion simulation) +. 0.5 *.
Mat.gemm_trace m_P (Mat.add m_Hc m_F)
On itere:
let rec iteration m_C n =
let m_P =
(* P = 2 C.C^t *)
gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C
in
let m_Hc, m_J, m_K =
let f =
Fock.make_rhf ~density:m_P ao_basis
in
Fock.(core f, coulomb f, exchange f)
in
let m_F =
Mat.add m_Hc (Mat.sub m_J m_K)
in
let m_C =
let m_C', _ =
Util.xt_o_x m_F m_X
|> Util.diagonalize_symm
in
gemm m_X m_C'
in
let energy =
(Simulation.nuclear_repulsion simulation) +. 0.5 *.
Mat.gemm_trace m_P (Mat.add m_Hc m_F)
in
Printf.printf "%f\n%!" energy;
if n > 0 then
iteration m_C (n-1)
let () = iteration m_C 20
Localisation des orbitales¶
Rotation des 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 $$
Construction de la matrice de rotation R de taille n par n:¶
Ainsi la construction d'une matrice de rotaion de taille n par n peut se faire de la manière suivante :
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^2 = 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.
Enfin on peut calculer R comme : $$ R = W \cos(\tau) W^\dagger + W \tau^{-1} \sin (\tau) W^\dagger X $$
(* Construction de la matrice de rotation R de taille n par n *)
Principe de la rotation par paires d'orbitales.¶
La rotation des orbitales se fait de la manière suivante :
On extrait de la matrice des OMs Phi 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 la matrice :
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 Psi (n par n) avec tous ses éléments nuls sauf les colonnes qui contiennent les OMs i et j
- une matrice Psi~ (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 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~.
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.
Localisation de Boyls¶
$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}2 = \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 (Boys)¶
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}
Localisation de Edmiston, C., & Ruedenberg, K.¶
(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)= \sum_n [Phi_n^2 | Phi_n^2 ]$
$= \sum_n < 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}^2 + B_{12}^2)^{1/2}$
$\alpha = (1/4) Acos (-A_{12} / (A_{12}^2 + B_{12}^2)^{1/2})$
$Sin(4 \alpha)= B_{12} / (A_{12}^2 + B_{12}^2)^{1/2}$
$\alpha = (1/4) Asin (B_{12} / (A_{12}^2 + B_{12}^2)^{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 ] $
Calculs : Localisation de Edminston¶
(* Localisation de Edminstion *)
let ee_ints = AOBasis.ee_ints ao_basis;;
let m_C = MOBasis.mo_coef mo_basis;;
let n_ao = Mat.dim1 m_C ;;
let sum a =
Array.fold_left (fun accu x -> accu +. x) 0. a
;;
(*Calcul des intégrales*)
let integral_general g i j =
Array.map (fun a ->
let v =
Array.map (fun b ->
let u =
Array.map (fun e ->
let t = Array.map (fun f ->
(g a b e f i j) *. ERI.get_phys ee_ints a e b f
) (Util.array_range 1 n_ao)
in sum t
) (Util.array_range 1 n_ao)
in sum u
) (Util.array_range 1 n_ao)
in sum v
) (Util.array_range 1 n_ao)
|> sum
(* Calcul de toutes les intégrales B_12 -> Matrice *)
let m_b12 = Mat.init_cols n_ao n_ao (fun i j ->
integral_general (fun a b e f i j ->
( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}
) i j
)
(* Calcul de toutes les intégrales A_12 -> Matrice *)
let m_a12 = Mat.init_cols n_ao n_ao (fun i j ->
integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,j}
-. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} )
*. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )
) i j
)
(* Calcul de tous les alpha -> Matrice *)
let m_alpha = Mat.init_cols n_ao n_ao ( fun i j ->
if i<>j then asin(m_b12.{i,j} /. sqrt((m_a12.{i,j})**2 + (m_b12.{i,j})**2 )
)
)
(* Extraction du plus grand angle alpha*)
let ij_max_alpha = Mat.as_vec m_alpha
|> iamax;;
(* Indices du plus grand angle alpha *)
let indice_i = ij_max_alpha mod n_ao;;
let indice_j = ij_max_alpha / n_ao + 1;;
Rotation des orbitales¶
let m_Phi = m_C
(* Matrice de rotation 2 par 2 *)
let m_R = 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 i et j de la matrice des OMs pour les mettres dans la matrice Ksi (n par 2)
pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)
let p = indice_i;;
let q = indice_j;;
let m_Ksi = Mat.init_cols n 2 (fun i j -> if j=1 then m_Phi.{i,p} else m_Phi.{i,q} );;
(* Fonction de calcul de ksi~ (matrice n par 2), nouvelle matrice par application de la matrice de rotation dans laquelle on obtient
les deux orbitales que l'on va réinjecter dans la matrice Phi*)
let m_Ksi_tilde = gemm ksi m_R;;
Réinjection¶
(* 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~ (i~ et j~) dans la matrice Phi *)
let m_Psi_tilde = Mat.init_cols n n (fun i j -> if j=q then m_Ksi_tilde.{i,p}
else if j=p then m_Ksi_tilde.{i,q}
else 0.);;
(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *)
let m_Psi = Mat.init_cols n n (fun i j -> if j=q then m_Ksi.{i,p}
else if j=p then m_Ksi.{i,q}
else 0.);;
(* 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
par la matrice *)
let m_interm = Mat.sub m_C M_Psi;;
(* Construction de la nouvelle matrice d'OMs phi~ *)
let m_Phi_tilde = Mat.add m_Psi_tilde m_interm;;
(* Fonction pour itérer *) (*cc -> critère de convergence*)
while m_alpha.{indice_i,indice_j} > cc do m_C = m_Phi_tilde
OLD ->¶
let ee_ints = AOBasis.ee_ints ao_basis;;
let m_C = MOBasis.mo_coef mo_basis;;
let n_ao = Mat.dim1 m_C ;;
let sum a =
Array.fold_left (fun accu x -> accu +. x) 0. a
;;
let i = 1;;
let j = 2;;
(*Calcul des intégrales*)
let integral_general g i j =
Array.map (fun a ->
let v =
Array.map (fun b ->
let u =
Array.map (fun e ->
let t = Array.map (fun f ->
(g a b e f i j) *. ERI.get_phys ee_ints a e b f
) (Util.array_range 1 n_ao)
in sum t
) (Util.array_range 1 n_ao)
in sum u
) (Util.array_range 1 n_ao)
in sum v
) (Util.array_range 1 n_ao)
|> sum
let func_a12 =
integral_general (fun a b e f i j ->
m_C.{a,i} *. m_C.{b,i} *. m_C.{e,i} *. m_C.{f,j}
-. 0.25 *. ( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} )
*. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} )
)
let func_b12 =
integral_general (fun a b e f i j ->
( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}
)
let i =1;;
let j = 2;;
let b12 =
integral_general (fun a b e f i j ->
( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}
)
(*Calcul de alpha*)
let alpha = asin(func_b12 /. sqrt((func_a12)**2. +. (func_b12)**2.));;
(*Matrice b12*)
let m_b12 = Mat.init_cols 10 10 (fun i j ->
integral_general (fun a b e f i j ->
( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}
) i j
)
let vec_test = Mat.as_vec m_b12
|> iamax;;
let test = Mat.max2 m_b12;;
$\langle \phi_i \phi_j | \phi_i \phi_j \rangle = \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} \sum_b \left(c_{bj} \sum_e \left(c_{ei} \sum_f \left(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) $
$\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 \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 ] $
$\langle \phi_i | x | \phi_j \rangle$ $=\langle \sum_a c_{ia} \chi_a | x | \sum_b c_{jb} \chi_b\rangle$ $=\sum_a \sum_b c_{jb} c_{ia} \langle \chi_a | x | \chi_b\rangle$
let ints = MOBasis.ee_ints mo_basis in
ERI.get_phys ints 1 2 1 2;;