2020-04-28 16:15:23 +02:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"source": [
"# Installation de QCaml"
]
},
{
"cell_type": "markdown",
"metadata": {
"hidden": true
},
"source": [
"1. Clonage de QCaml:\n",
" ```bash\n",
" git clone https://gitlab.com/scemama/QCaml.git\n",
" ```\n",
"2. Installation des dependances:\n",
" ```bash\n",
" opam install ocamlbuild ocamlfind lacaml gnuplot getopt alcotest zarith\n",
" cd QCaml\n",
" ./configure\n",
" ```\n",
"3. Compilation\n",
" ```bash\n",
" make\n",
" ```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Liens utiles"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Documentation des modules : https://sanette.github.io/ocaml-api/4.10/index.html\n",
"* Lacaml : http://mmottl.github.io/lacaml/api/lacaml/\n",
"* Gnuplot : https://github.com/c-cube/ocaml-gnuplot"
]
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"source": [
"# Initialisation"
]
},
{
"cell_type": "markdown",
"metadata": {
"hidden": true
},
"source": [
"Bloc a executer avant de pouvoir utiliser QCaml dans le Notebook"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let png_image = print_endline ;;\n",
"\n",
"(* --------- *)\n",
"\n",
"(*Mettre le bon chemin ici *)\n",
2020-05-14 16:06:31 +02:00
"#cd \"/home/ydamour/QCaml\";;\n",
2020-04-28 16:15:23 +02:00
"\n",
"#use \"topfind\";;\n",
"#require \"jupyter.notebook\";;\n",
2020-05-01 08:45:25 +02:00
" \n",
2020-04-28 16:15:23 +02:00
"#require \"gnuplot\";;\n",
"let png_image name = \n",
" Jupyter_notebook.display_file ~base64:true \"image/png\" (\"Notebooks/images/\"^name)\n",
";;\n",
"\n",
"#require \"lacaml.top\";;\n",
"#require \"alcotest\";;\n",
"#require \"str\";;\n",
"#require \"bigarray\";;\n",
"#require \"zarith\";;\n",
"#require \"getopt\";;\n",
"#directory \"_build\";;\n",
"#directory \"_build/Basis\";;\n",
"#directory \"_build/CI\";;\n",
"#directory \"_build/MOBasis\";;\n",
"#directory \"_build/Nuclei\";;\n",
"#directory \"_build/Parallel\";;\n",
"#directory \"_build/Perturbation\";;\n",
"#directory \"_build/SCF\";;\n",
"#directory \"_build/Utils\";;\n",
"\n",
"\n",
"#load \"Constants.cmo\";;\n",
"#load_rec \"Util.cma\";;\n",
"#load_rec \"Matrix.cmo\";;\n",
"#load_rec \"Simulation.cmo\";;\n",
"#load_rec \"Spindeterminant.cmo\";;\n",
"#load_rec \"Determinant.cmo\";;\n",
"#load_rec \"HartreeFock.cmo\";;\n",
"#load_rec \"MOBasis.cmo\";;\n",
"#load_rec \"F12CI.cmo\";;\n",
"\n",
"#install_printer AngularMomentum.pp_string ;;\n",
"#install_printer Basis.pp ;;\n",
"#install_printer Charge.pp ;;\n",
"#install_printer Coordinate.pp ;;\n",
"#install_printer Vector.pp;;\n",
"#install_printer Matrix.pp;;\n",
"#install_printer Util.pp_float_2darray;;\n",
"#install_printer Util.pp_float_array;;\n",
"#install_printer Util.pp_matrix;;\n",
"#install_printer HartreeFock.pp ;;\n",
"#install_printer Fock.pp ;;\n",
"#install_printer MOClass.pp ;;\n",
"#install_printer DeterminantSpace.pp;;\n",
"#install_printer SpindeterminantSpace.pp;;\n",
"#install_printer Bitstring.pp;;\n",
"let pp_mo ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;\n",
"#install_printer pp_mo;;\n",
"\n",
"\n",
"(* --------- *)\n",
"\n",
"open Lacaml.D\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Calculs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## H$_2$O"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
2020-05-15 12:33:33 +02:00
"(*\n",
2020-04-28 16:15:23 +02:00
"let xyz_string = \n",
"\"3\n",
"Water\n",
"O 0. 0. 0.\n",
"H -0.756950272703377558 0. -0.585882234512562827\n",
"H 0.756950272703377558 0. -0.585882234512562827\n",
2020-05-15 12:33:33 +02:00
"\"\n",
"*)"
2020-05-01 08:45:25 +02:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
2020-05-07 21:07:04 +02:00
"### H$_4$"
2020-05-01 08:45:25 +02:00
]
},
2020-05-07 21:07:04 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-08 10:49:21 +02:00
"(* Fonction création chaîne linéaire de n H *)\n",
2020-05-15 12:33:33 +02:00
"\n",
2020-05-08 10:49:21 +02:00
"let xyz d n = \n",
2020-05-07 21:07:04 +02:00
" let accu = \"\"\n",
" in\n",
2020-05-08 10:49:21 +02:00
" let rec toto accu d n =\n",
2020-05-07 21:07:04 +02:00
" let accu = \n",
" if n=0 \n",
" then accu ^ \"\"\n",
" else\n",
" match n with \n",
2020-05-08 10:49:21 +02:00
" | 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",
2020-05-07 21:07:04 +02:00
" in\n",
" accu\n",
2020-05-08 10:49:21 +02:00
" in string_of_int(n) ^ \"\\nH\" ^ string_of_int(n) ^ \"\\n\" ^ toto accu d n;;\n",
2020-05-07 21:07:04 +02:00
" \n",
2020-05-15 16:31:00 +02:00
"let xyz_string = xyz 1.8 6;;\n"
2020-04-28 16:15:23 +02:00
]
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"source": [
"### Base atomique"
]
},
{
"cell_type": "markdown",
"metadata": {
"hidden": true
},
"source": [
"Les bases atomiques sont telechargeables ici: https://www.basissetexchange.org/\n",
"\n",
"On telecharge la base cc-pvdz depuis le site:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
2020-05-15 12:33:33 +02:00
"\n",
2020-05-14 08:54:39 +02:00
"let basis_string = \n",
2020-05-07 18:00:34 +02:00
"\"\n",
2020-05-01 08:45:25 +02:00
"HYDROGEN\n",
2020-05-07 18:00:34 +02:00
"S 6\n",
"1 0.3552322122E+02 0.9163596281E-02\n",
"2 0.6513143725E+01 0.4936149294E-01\n",
"3 0.1822142904E+01 0.1685383049E+00\n",
"4 0.6259552659E+00 0.3705627997E+00\n",
"5 0.2430767471E+00 0.4164915298E+00\n",
"6 0.1001124280E+00 0.1303340841E+00\n",
2020-05-01 08:45:25 +02:00
"\"\n",
2020-05-07 18:00:34 +02:00
"\n",
"\n",
"\n",
2020-05-15 12:33:33 +02:00
"(*\n",
2020-05-13 15:46:48 +02:00
"let basis_string = \n",
"\"\n",
2020-04-28 16:15:23 +02:00
"HYDROGEN\n",
"S 4\n",
"1 1.301000E+01 1.968500E-02\n",
"2 1.962000E+00 1.379770E-01\n",
"3 4.446000E-01 4.781480E-01\n",
"4 1.220000E-01 5.012400E-01\n",
"S 1\n",
"1 1.220000E-01 1.000000E+00\n",
"P 1\n",
"1 7.270000E-01 1.0000000\n",
2020-05-14 08:54:39 +02:00
"\n",
2020-04-28 16:15:23 +02:00
"OXYGEN\n",
"S 9\n",
"1 1.172000E+04 7.100000E-04\n",
"2 1.759000E+03 5.470000E-03\n",
"3 4.008000E+02 2.783700E-02\n",
"4 1.137000E+02 1.048000E-01\n",
"5 3.703000E+01 2.830620E-01\n",
"6 1.327000E+01 4.487190E-01\n",
"7 5.025000E+00 2.709520E-01\n",
"8 1.013000E+00 1.545800E-02\n",
"9 3.023000E-01 -2.585000E-03\n",
"S 9\n",
"1 1.172000E+04 -1.600000E-04\n",
"2 1.759000E+03 -1.263000E-03\n",
"3 4.008000E+02 -6.267000E-03\n",
"4 1.137000E+02 -2.571600E-02\n",
"5 3.703000E+01 -7.092400E-02\n",
"6 1.327000E+01 -1.654110E-01\n",
"7 5.025000E+00 -1.169550E-01\n",
"8 1.013000E+00 5.573680E-01\n",
"9 3.023000E-01 5.727590E-01\n",
"S 1\n",
"1 3.023000E-01 1.000000E+00\n",
"P 4\n",
"1 1.770000E+01 4.301800E-02\n",
"2 3.854000E+00 2.289130E-01\n",
"3 1.046000E+00 5.087280E-01\n",
"4 2.753000E-01 4.605310E-01\n",
"P 1\n",
"1 2.753000E-01 1.000000E+00\n",
"D 1\n",
"1 1.185000E+00 1.0000000\n",
2020-05-15 12:33:33 +02:00
"\"\n",
"*)"
2020-04-28 16:15:23 +02:00
]
},
{
"cell_type": "markdown",
"metadata": {
"hidden": true
},
"source": [
"Une orbitale atomique centree sur l'atome A est composee d'une contraction de m Gaussiennes:\n",
"$$\n",
"\\chi(r) = \\sum_{i=1}^m a_i \\exp \\left( -\\alpha_i |r-r_A|^2 \\right)\n",
"$$\n",
"Dans le fichier de base, la 2ieme colonne represente l'exposant $\\alpha_i$ et la 3ieme colonne le coefficient de contraction $a_i$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let nuclei =\n",
" Nuclei.of_xyz_string xyz_string\n",
" \n",
"let basis = \n",
" Basis.of_nuclei_and_basis_string nuclei basis_string\n",
" \n",
"let simulation = \n",
" Simulation.make ~charge:0 ~multiplicity:1 ~nuclei basis\n",
" \n",
"let ao_basis = \n",
2020-05-15 17:14:02 +02:00
" Simulation.ao_basis simulation\n",
" \n",
"let nocc =\n",
" let elec = Simulation.electrons simulation in\n",
" Electrons.n_alfa elec\n"
2020-04-28 16:15:23 +02:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot des orbitales atomiques"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let plot data filename = \n",
" let output = Gnuplot.Output.create (`Png filename) in\n",
" let gp = Gnuplot.create () in\n",
" Gnuplot.set gp ~use_grid:true;\n",
" List.map Gnuplot.Series.lines_xy data\n",
" |> Gnuplot.plot_many gp ~output;\n",
" Gnuplot.close gp ;\n",
" Jupyter_notebook.display_file ~base64:true \"image/png\" filename\n",
";;\n",
"\n",
"let x_values = \n",
" let n = 1000 in\n",
" \n",
" let xmin, xmax =\n",
" let coord =\n",
" Array.map snd nuclei\n",
" |> Array.map (fun a -> Coordinate.(get X) a)\n",
" in\n",
" Array.sort compare coord;\n",
" coord.(0) -. 4. ,\n",
" coord.(Array.length coord -1) +. 4.\n",
" in\n",
"\n",
" let dx =\n",
" (xmax -. xmin) /. (float_of_int n -. 1.)\n",
" in\n",
" Array.init n (fun i -> xmin +. (float_of_int i)*.dx)\n",
"in\n",
"\n",
"let data = \n",
" Array.map (fun x -> \n",
" let point = Coordinate.make_angstrom\n",
" { Coordinate.\n",
" x ; y = 0. ; z = 0.\n",
" } in\n",
" AOBasis.values ao_basis point\n",
" ) x_values\n",
" |> Mat.of_col_vecs\n",
" |> Mat.transpose_copy\n",
" |> Mat.to_col_vecs\n",
" |> Array.map Vec.to_list\n",
" |> Array.map (fun l -> List.mapi (fun i y -> (x_values.(i),y)) l)\n",
" |> Array.to_list\n",
"in\n",
"plot data \"test_data.png\""
]
},
{
"cell_type": "markdown",
"metadata": {
"heading_collapsed": true
},
"source": [
"### Calcul Hartree-Fock"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let hf = HartreeFock.make simulation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let mo_basis = MOBasis.of_hartree_fock hf"
]
},
{
"cell_type": "markdown",
"metadata": {
"hidden": true
},
"source": [
"Orbitales moleculaires :\n",
"$$\n",
"\\phi_j(r) = \\sum_{i=1}^{N_b} C_{ij} \\chi_i(r)\n",
"$$\n",
"\n",
"* $i$: lignes\n",
"* $j$: colonnes"
]
},
{
"cell_type": "markdown",
"metadata": {
"hidden": true
},
"source": [
"Extraction des OM de la structure de donnees `mo_basis` comme une matrice $C$ utilisable avec Lacaml:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"hidden": true
},
"outputs": [],
"source": [
"let mo_coef = MOBasis.mo_coef mo_basis"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let plot data filename = \n",
" let output = Gnuplot.Output.create (`Png filename) in\n",
" let gp = Gnuplot.create () in\n",
" Gnuplot.set gp ~use_grid:true;\n",
" List.map Gnuplot.Series.lines_xy data\n",
" |> Gnuplot.plot_many gp ~output;\n",
" Gnuplot.close gp ;\n",
" Jupyter_notebook.display_file ~base64:true \"image/png\" filename\n",
";;\n",
"\n",
"let x_values = \n",
" let n = 1000 in\n",
" \n",
" let xmin, xmax =\n",
" let coord =\n",
" Array.map snd nuclei\n",
" |> Array.map (fun a -> Coordinate.(get X) a)\n",
" in\n",
" Array.sort compare coord;\n",
" coord.(0) -. 4. ,\n",
" coord.(Array.length coord -1) +. 4.\n",
" in\n",
"\n",
" let dx =\n",
" (xmax -. xmin) /. (float_of_int n -. 1.)\n",
" in\n",
" Array.init n (fun i -> xmin +. (float_of_int i)*.dx)\n",
"in\n",
"\n",
"let data = \n",
" let result = \n",
" Array.map (fun x -> \n",
" let point = Coordinate.make_angstrom\n",
" { Coordinate.\n",
" x ; y = 0. ; z = 0.\n",
" } in\n",
" MOBasis.values mo_basis point\n",
" ) x_values\n",
" |> Mat.of_col_vecs\n",
" |> Mat.transpose_copy\n",
" |> Mat.to_col_vecs\n",
" |> Array.map Vec.to_list\n",
" |> Array.map (fun l -> List.mapi (fun i y -> (x_values.(i),y)) l)\n",
" in\n",
" [ result.(0) ; result.(1) ; result.(2) ]\n",
"in\n",
"\n",
"plot data \"test_data.png\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calcul Hartree-Fock a la main"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Methode\n",
"\n",
"$$\n",
"\\hat{H} = \\hat{T} + \\hat{V}^{\\text{NN}} + \\hat{V}^{\\text{eN}} + \\hat{V}^{\\text{ee}}\n",
"$$\n",
"\n",
"On exprime les differentes quantites dans la base des orbitales atomiques pour chacun des termes:\n",
"\n",
"\\begin{eqnarray}\n",
"T_{ij} & = & \\langle \\chi_i | \\hat{T} | \\chi_j \\rangle =\n",
" \\iiint \\chi_i(\\mathbf{r}) \\left( -\\frac{1}{2} \\Delta\n",
" \\chi_j(\\mathbf{r}) \\right) d\\mathbf{r} \\\\\n",
"V^{NN} & = & \\text{constante} \\\\\n",
"V^{eN}_{ij} & = & \\langle \\chi_i | \\hat{V}^{\\text{eN}} | \\chi_j \\rangle =\n",
" \\sum_A \\iiint \\chi_i(\\mathbf{r}) \n",
" \\frac{-Z_A}{|\\mathbf{r} - \\mathbf{R}_A|}\n",
" \\chi_j(\\mathbf{r}) d\\mathbf{r} \\\\\n",
"V^{ee}_{ijkl} & = & \\langle \\chi_i \\chi_j | \\hat{V}^{\\text{ee}} | \\chi_k \\chi_l \\rangle =\n",
" \\iiint \\iiint \\chi_i(\\mathbf{r}_1) \\chi_j(\\mathbf{r}_2) \n",
" \\frac{1}{|\\mathbf{r}_1 - \\mathbf{r}_2|}\n",
" \\chi_k(\\mathbf{r}) \\chi_l(\\mathbf{r}) d\\mathbf{r}_1\n",
" d\\mathbf{r}_2 \\\\\n",
"\\end{eqnarray}\n",
"ou $\\mathbf{R}_A$ est la position du noyau $A$ et $Z_A$ sa charge.\n",
"\n",
"\n",
"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.\n",
"\n",
"On definit donc 2 operateurs pour chaque orbitale:\n",
"- Coulomb :\n",
"$$ \\hat{J}_j \\chi_i(\\mathbf{r}_1) = \\chi_i(\\mathbf{r}_1) \n",
" \\iiint \\frac{1}{|\\mathbf{r}_1 - \\mathbf{r}_2|} |\\chi_j(\\mathbf{r}_2)|^2 d \\mathbf{r}_2\n",
"$$\n",
"- Echange :\n",
"$$ \\hat{K}_j \\chi_i(\\mathbf{r}_1) = \\chi_j(\\mathbf{r}_1) \n",
" \\iiint \\frac{1}{|\\mathbf{r}_1 - \\mathbf{r}_2|} \\chi_i(\\mathbf{r}_2) \\chi_j(\\mathbf{r}_2) d \\mathbf{r}_2\n",
"$$\n",
"et on n'a plus que des operateurs a 1 electron.\n",
"\n",
"\\begin{eqnarray}\n",
"J_{ij} & = & \\sum_{kl} P_{kl} \\langle i k | j l \\rangle \\\\\n",
"K_{il} & = & \\sum_{kj} P_{kj} \\langle i k | j l \\rangle \n",
"\\end{eqnarray}\n",
"ou $P$ est la matrice densite definie comme\n",
"$$\n",
"P_{ij} = \\sum_{k=1}^{N_{\\text{occ}}} 2 C_{ik} C_{kj}\n",
"$$\n",
"et $C$ est la matrice des coefficients des orbitales moleculaires exprimees dans la base des orbitales atomiques.\n",
"\n",
"Une orbitale moleculaire est une combinaison lineaire d'orbitale atomique:\n",
"$$\n",
"\\phi_k(\\mathbf{r}) = \\sum_i C_{ik} \\chi_i(\\mathbf{r})\n",
"$$\n",
"\n",
"\n",
"Les orbitales moleculaires sont toutes orthogonales entre elles, et normalisees a 1.\n",
"\n",
"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:\n",
"$$\n",
"\\Psi(\\mathbf{r}_1,\\dots,\\mathbf{r}_N) = \n",
"\\left| \\begin{array}{ccc}\n",
"\\phi_1(\\mathbf{r}_1) & \\dots & \\phi_N(\\mathbf{r}_1) \\\\\n",
"\\vdots & \\ddots & \\vdots \\\\\n",
"\\phi_1(\\mathbf{r}_N) & \\dots & \\phi_N(\\mathbf{r}_N) \n",
"\\end{array}\\right| \n",
"$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Application"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On commence par creer une base orthogonale $X$ a partir des OA"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let m_X = \n",
" AOBasis.ortho ao_basis\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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}$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let c_of_h m_H = \n",
" (* On exprime H dans la base orthonormale *)\n",
" let m_Hmo =\n",
" Util.xt_o_x m_H m_X (* H_mo = X^t H X *)\n",
" in\n",
" \n",
" (* On diagonalise cet Hamiltonien *)\n",
" let m_C', _ =\n",
" Util.diagonalize_symm m_Hmo\n",
" in\n",
" \n",
" (* On re-exprime les MOs dans la base des AOs (non-orthonormales) *) \n",
" gemm m_X m_C' (* C = X.C' *)\n",
" \n",
" \n",
"let m_C = \n",
" match Guess.make ~nocc ~guess:`Hcore ao_basis with\n",
" | Hcore m_H -> c_of_h m_H\n",
" | _ -> assert false\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"--------\n",
"On construit la matrice densite"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let m_P = \n",
" (* P = 2 C.C^t *)\n",
" gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On construit les matrices $H_c = T+V_{\\text{eN}}$, et $J$ et $K$:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let m_Hc, m_J, m_K =\n",
" let f =\n",
" Fock.make_rhf ~density:m_P ao_basis\n",
" in\n",
" Fock.(core f, coulomb f, exchange f)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On construit l'operateur de Fock:\n",
"$ F = H_c + J - K $"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let m_F = \n",
" Mat.add m_Hc (Mat.sub m_J m_K)\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On l'exprime dans la base orthonormale, on le diagonalise, et on le re-exprime ses vecteurs propres dans la base d'OA."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let m_C =\n",
" let m_C', _ = \n",
" Util.xt_o_x m_F m_X\n",
" |> Util.diagonalize_symm\n",
" in\n",
" gemm m_X m_C'\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On calcule l'energie avec ces nouvelles orbitales:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let energy =\n",
" (Simulation.nuclear_repulsion simulation) +. 0.5 *.\n",
" Mat.gemm_trace m_P (Mat.add m_Hc m_F)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On itere:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let rec iteration m_C n =\n",
" let m_P = \n",
" (* P = 2 C.C^t *)\n",
" gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C\n",
" in\n",
"\n",
" let m_Hc, m_J, m_K =\n",
" let f =\n",
" Fock.make_rhf ~density:m_P ao_basis\n",
" in\n",
" Fock.(core f, coulomb f, exchange f)\n",
" in\n",
"\n",
" let m_F = \n",
" Mat.add m_Hc (Mat.sub m_J m_K)\n",
" in\n",
"\n",
" let m_C =\n",
" let m_C', _ = \n",
" Util.xt_o_x m_F m_X\n",
" |> Util.diagonalize_symm\n",
" in\n",
" gemm m_X m_C'\n",
" in\n",
"\n",
" let energy =\n",
" (Simulation.nuclear_repulsion simulation) +. 0.5 *.\n",
" Mat.gemm_trace m_P (Mat.add m_Hc m_F)\n",
" in\n",
" Printf.printf \"%f\\n%!\" energy;\n",
" if n > 0 then\n",
" iteration m_C (n-1)\n",
"\n",
"let () = iteration m_C 20\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Localisation des orbitales"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Rotation des orbitales\n",
"\n",
"On peut écrire une matrice de rotation $R$ comme\n",
"$$\n",
"R = \\exp(X)\n",
"$$\n",
"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:\n",
"$$\n",
"R = \\left( \\begin{array}{cc}\n",
"\\cos \\theta & -\\sin \\theta \\\\\n",
"\\sin \\theta & \\cos \\theta \\end{array} \\right)\n",
"= \\exp\n",
" \\left( \\begin{array}{cc}\n",
" 0 & \\theta \\\\\n",
"- \\theta & 0 \\end{array} \\right)\n",
"$$\n",
"\n",
"D'abord on diagonalise $X^2$ :\n",
"$$\n",
"X^2 = W(-\\tau^2) W^\\dagger.\n",
"$$\n",
"Ensuite, $R$ peut être écrit comme\n",
"$$\n",
"R = W \\cos(\\tau) W^\\dagger + W \\tau^{-1} \\sin (\\tau) W^\\dagger X\n",
"$$\n",
"\n",
"## Construction de la matrice de rotation R de taille n par n:\n",
"\n",
"Ainsi la construction d'une matrice de rotaion de taille n par n peut se faire de la manière suivante :\n",
"\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).\n",
"\n",
"Tout d'abord on met au carré cette matrice X ( gemm X X ) tel que : \n",
" \n",
"$$\n",
"X^2 = X X\n",
"$$\n",
" \n",
"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.\n",
"\n",
"Enfin on peut calculer R comme :\n",
"$$\n",
"R = W \\cos(\\tau) W^\\dagger + W \\tau^{-1} \\sin (\\tau) W^\\dagger X\n",
"$$\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Construction de la matrice de rotation R de taille n par n *)\n"
]
},
{
"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 OMs Phi 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 la matrice :\n",
"\n",
"R = ( cos(alpha) -sin(alpha) )\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~. 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 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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
2020-05-04 16:14:06 +02:00
"## Localisation de Boys"
2020-04-28 16:15:23 +02:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$N$ orthonormal molecular orbitals\n",
"$$\n",
"\\phi_i({\\bf r})= \\sum_{k=1}^N c_{ik} \\chi_k\n",
"$$\n",
"$$\n",
"{\\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\n",
"$$\n",
"$$\n",
"{\\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]\n",
"$$\n",
"$$\n",
"{\\cal B}_2 = \\sum_{i=1}^N \\big[ \\langle x^4 \\rangle_i - 4 \\langle x^3 \\rangle_i \\langle x \\rangle_i\n",
" + 6 \\langle x^2 \\rangle_i \\langle x \\rangle^2_i\n",
"- 3 \\langle x \\rangle^4_i \\big] + \\big[ ...y...] + \\big[ ...z...] \n",
"$$\n",
"Minimization of ${\\cal B}$ with respect to an arbitrary rotation $R$\n",
"$$\n",
" \\langle R \\phi_i x^n R \\phi_i \\rangle = \\sum_{k,l=1}^N R_{ik} R_{il} \\langle \\phi_k| x^n | \\phi_l …\\rangle= \n",
"$$\n",
"$$\n",
"\\sum_{k,l=1}^N R_{ik} R_{il} \\sum_{m,o=1}^N c_{km} c_{ln} \\langle \\chi_m | x^n |\\chi_o \\rangle\n",
"$$\n",
"We need to compute\n",
"$$\n",
"S^x_{n;mo}= \\langle \\chi_m | x^n |\\chi_o \\rangle\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Rotation de deux orbitales (Boys)"
]
},
2020-04-29 17:21:24 +02:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On part de la matrice des OMs donné par $HF$ -> m_C.On calcul $A_{ij}^{x,y,z}$, et $B_{ij}^{x,y,z}$, tel que :\n",
"\n",
"$A_{ij}^x = \\langle \\phi_i | x^2 | \\phi_j \\rangle$ \n",
"\n",
"$= \\sum_a \\sum_b c_{ai} c_{bj} \\langle \\chi_a | x^2 | \\chi_b \\rangle$\n",
"\n",
"$B_{ij}^x = \\langle \\phi_i | x | \\phi_j \\rangle$ \n",
"\n",
"$= \\sum_a \\sum_b c_{ai} c_{bj} \\langle \\chi_a | x | \\chi_b \\rangle$\n",
"\n",
"On forme ainsi des matrices (n par n) $A_{ij}^x$, $A_{ij}^y$, $A_{ij}^z$, $B_{ij}^x$, $B_{ij}^y$ et $B_{ij}^z$ que l'on peut sommer pour obtenir :\n",
"\n",
"$A_{ij} =A_{ij}^x + A_{ij}^y + A_{ij}^z$ et $B_{ij} =B_{ij}^x + B_{ij}^y + B_{ij}^z$\n",
"\n",
"$D$ est défini comme :\n",
"\n",
"$D(\\theta) = D(0) + \\frac{1}{4} [(1-cos4\\theta)\\beta_{ij}+sin4\\theta \\gamma_{ij}] $\n",
"\n",
"Où\n",
"\n",
"$\\beta_{ij}= (B^x_{ii}-B^x_{jj})^2 - 4 {(B^x_{ij})}^2 + [...y...] + [...z...]$\n",
"\n",
"$=(B_{ii}-B_{jj})^2 - 4 {(B_{ij})}^2 $\n",
"\n",
"Et \n",
"\n",
"$\\gamma_{ij}= 4 B^x_{ij} (B^{x}_{ii}-B^x_{jj}) + [...y...] + [...z...]$\n",
"\n",
"$=4 B_{ij} (B_{ii}-B_{jj})$\n",
"\n",
"Avec\n",
"\n",
"$D(0)= D^{x}(0) + D^{y}(0) + D^{z}(0)$\n",
"\n",
"$=A_{ii}^x + A_{jj}^x - (\\tilde B_{ii}^x)^2 - (\\tilde B_{jj}^x)^2$\n",
"$+A_{ii}^y + A_{jj}^y - (\\tilde B_{ii}^y)^2 - (\\tilde B_{jj}^y)^2$\n",
"$+A_{ii}^z + A_{jj}^z - (\\tilde B_{ii}^z)^2 - (\\tilde B_{jj}^z)^2$\n",
"\n",
"$= A_{ii} + A_{jj} - \\tilde B_{ii}^2 - \\tilde B_{jj}^2$\n",
"\n",
"Avec\n",
"\n",
"${\\tilde B}^{x2}_{ii}= (cos^2\\theta B^x_{ii} + sin^2\\theta B^x_{jj} - sin2\\theta B^x_{ij})^2$\n",
"\n",
"${\\tilde B}^{x2}_{jj}= (sin^2\\theta B^x_{ii} + cos^2\\theta B^x_{jj} + sin2\\theta B^x_{ij})^2$\n",
"\n",
"Et donc :\n",
"\n",
"${\\tilde B}^{2}_{ii}= (cos^2\\theta B_{ii} + sin^2\\theta B_{jj} - sin2\\theta B_{ij})^2$\n",
"\n",
"${\\tilde B}^{2}_{jj}= (sin^2\\theta B_{ii} + cos^2\\theta B_{jj} + sin2\\theta B_{ij})^2$\n",
"\n",
"Ainsi, on a : \n",
"\n",
"$D(\\theta) = A_{ii} + A_{jj} - (cos^2\\theta B_{ii} + sin^2\\theta B_{jj} - sin2\\theta B_{ij})^2 - (sin^2\\theta B_{ii} + cos^2\\theta B_{jj} + sin2\\theta B_{ij})^2$\n",
"\n",
2020-05-01 08:45:25 +02:00
"Les extrema de $D(\\theta)$ sont de la forme :\n",
2020-04-30 11:07:37 +02:00
"\n",
"$Tan(\\theta) = - \\frac{\\gamma}{\\beta}$\n",
"\n",
2020-05-01 08:45:25 +02:00
"$D(\\theta)$ à 4 extrema :\n",
2020-04-30 11:07:37 +02:00
"\n",
"$4\\theta; \\;\\; 4\\theta +\\pi; \\;\\; 4\\theta+ 2\\pi; \\;\\; 4\\theta+ 3\\pi$\n",
"\n",
2020-05-01 08:45:25 +02:00
"$\\theta_{ij}= Atan(- \\frac{\\gamma}{\\beta})$\n",
2020-04-29 17:21:24 +02:00
"\n",
2020-05-01 08:45:25 +02:00
"Ainsi on peut calculer la matrice des $\\theta$ et effectuer la rotation pour le couple d'orbitale $(i,j)$ ayant le $\\theta_{max}$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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",
"$- \\frac {1}{4}(\\langle \\phi_1 | \\bar r | \\phi_1 \\rangle $\n",
2020-05-16 11:09:53 +02:00
"$- \\langle \\phi_2 | \\bar r | \\phi_2 \\rangle . \\langle \\phi_1 | \\bar r | \\phi_1 \\rangle$\n",
2020-05-01 08:45:25 +02:00
"$- \\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",
2020-05-16 11:09:53 +02:00
"$ . \\langle \\phi_1 | \\bar r | \\phi_2 \\rangle $\n",
2020-05-01 08:45:25 +02:00
"\n",
"Avec \n",
"\n",
"$A^r_{12}=A^x_{12} + A^y_{12} + A^z_{12}$\n",
"\n",
2020-05-08 10:49:21 +02:00
"$B^r_{12}=B^x_{12} + B^y_{12} + B^z_{12}$\n",
"\n",
2020-05-16 11:09:53 +02:00
"Et le critère à maximiser est :\n",
2020-05-08 10:49:21 +02:00
"\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"
2020-04-29 17:21:24 +02:00
]
},
2020-04-28 16:15:23 +02:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Rotation of angle $\\theta$\n",
"$$\n",
"\\tilde{\\phi}_1 = cos\\theta \\phi_1 -sin\\theta \\phi_2\n",
"$$\n",
"$$\n",
"\\tilde{\\phi}_2 = sin\\theta \\phi_1 +cos\\theta \\phi_2\n",
"$$\n",
"Let us note\n",
"$$\n",
"{\\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}}\n",
"$$\n",
"and\n",
"\\begin{eqnarray}\n",
"A^x_{ij} & = & \\langle \\phi_i| x^2 | \\phi_j \\rangle \\\\\n",
"B^x_{ij} & = & \\langle \\phi_i| x | \\phi_j \\rangle \\;\\; i,j=1,2 \n",
"\\end{eqnarray}\n",
"We have\n",
"$$\n",
"{\\cal B}_x(\\theta) = A^x_{11}+A^x_{22} - ({\\tilde B}^{x2}_{11} + {\\tilde B}^{x2}_{22})\n",
"$$\n",
"with\n",
"$$\n",
"{\\tilde B}^{x2}_{11}= (cos^2\\theta B^x_{11} + sin^2\\theta B^x_{22} - sin2\\theta B^x_{12})^2\n",
"$$\n",
"$$\n",
"{\\tilde B}^{x2}_{22}= (sin^2\\theta B^x_{11} + cos^2\\theta B^x_{22} + sin2\\theta B^x_{12})^2\n",
"$$\n",
"$$\n",
"{\\cal B}(\\theta)= {\\cal B}_x(\\theta)+{\\cal B}_y(\\theta)+{\\cal B}_z(\\theta)=(x) + (y)+(z)\n",
"$$\n",
"with\n",
"$$\n",
"{\\cal B}_x(\\theta)= A^x_{11}+A^x_{22} \n",
"- [(cos^4\\theta+ sin^4\\theta)({B^x_{11}}^2+ {B^x_{22}}^2 )\n",
"+ 2 sin^2 2\\theta {B^x_{12}}^2\n",
"+ 2 sin 2\\theta cos 2\\theta (({B^x_{22}} -{B^x_{11}} ) {B^x_{12}}]\n",
"$$\n",
"and idem for $y$ and $z$. \n",
"Using the fact that\n",
"$$\n",
"cos^4\\theta+ sin^4\\theta= \\frac{1}{4} ( 3 + cos4\\theta)\n",
"$$\n",
"$$\n",
"{\\cal B}_x(\\theta)= A^x_{11}+A^x_{22} \n",
"- [ \\frac{1}{4} ( 3 + cos4\\theta)({B^x_{11}}^2+ {B^x_{22}}^2 )\n",
"+ (1 -cos 4\\theta) {B^x_{12}}^2\n",
"+ sin 4\\theta (({B^x_{22}} -{B^x_{11}} ) {B^x_{12}}]\n",
"$$\n",
"Finally, we get\n",
"\\begin{equation}\n",
"{\\cal B}(\\theta)= {\\cal B}(0) + \\frac{1}{4} [(1-cos4\\theta)\\beta+sin4\\theta \\gamma] \n",
"\\end{equation}\n",
"where\n",
"$$\n",
"{\\cal B}(0)= A^x_{11}+A^x_{22} -((B^{x}_{11})^2+(B^{x}_{22})^2) + [...y...] + [...z...]\n",
"$$\n",
"$$\n",
"\\beta= (B^x_{11}-B^x_{22})^2 - 4 {(B^x_{12})}^2 + [...y...] + [...z...]\n",
"$$\n",
"and\n",
"$$\n",
"\\gamma= 4 B^x_{12} (B^{x}_{11}-B^x_{22}) + [...y...] + [...z...]\n",
"$$\n",
"Let us compute the derivative; we get\n",
"$$\n",
"\\frac{\\partial {\\cal B}(\\theta)}{\\partial \\theta} = \n",
"\\beta sin4\\theta \n",
"+ \\gamma cos4\\theta\n",
"$$\n",
"Extrema of ${\\cal B}(\\theta)$\n",
"\\begin{equation}\n",
"tg4\\theta= -\\frac{\\gamma}{\\beta} \n",
"\\end{equation}\n",
"There are four extrema:\n",
"$$\n",
"4\\theta; \\;\\; 4\\theta +\\pi; \\;\\; 4\\theta+ 2\\pi; \\;\\; 4\\theta+ 3\\pi\n",
"$$\n",
"Value of the second derivative of $\\cal{B}$ at the extrema\n",
"Value of $\\cal{B}$ at the extrema\n",
"\\begin{equation}\n",
"\\frac{\\partial^2 B(\\theta)}{\\partial \\theta^2}= 4 cos4\\theta \\frac{\\beta^2 + \\gamma^2}{\\beta}\n",
"\\end{equation}\n",
"There are two minima and two maxima since $cos4\\theta= -cos4(\\theta+\\pi)= -cos4(\\theta+2\\pi)=-cos4(\\theta+3\\pi)$.\n",
"Value of $\\cal{B}$ at the extrema\n",
"\\begin{equation}\n",
"{\\cal B}(\\theta)= {\\cal B}(0) + \\frac{1}{4} (\\beta -\\frac{\\beta^2 + \\gamma^2}{\\beta} {cos4\\theta})\n",
"\\end{equation}"
]
},
{
"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–
"\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²_n | 1/r_{12} | Phi²_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 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",
"\n",
"$j~ (x) = -sin(\\gamma) i(x) + cos(\\gamma) j(x)$\n",
"\n",
2020-05-06 15:39:23 +02:00
"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",
2020-04-28 16:15:23 +02:00
"\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",
2020-05-08 10:49:21 +02:00
"\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",
2020-04-28 16:15:23 +02:00
" \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",
2020-04-29 17:21:24 +02:00
"$= \\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",
"On extrait de la matrice des OMs les orbitales i et j (vecteurs colonnes) pour former une nouvelle matrice m_Ksi de dimension n * 2. \n",
"\n",
"Pour effectuer la rotation des orbitales i et j on utilise la matrice de rotation m_R pour la rotation de 2 orbitales, qui est définie comme :\n",
"\n",
"m_R =\n",
"\n",
" ( cos(alpha) -sin(alpha) )\n",
"\n",
" ( sin(alpha) cos(alpha) )\n",
" \n",
"On applique m_R à m_Ksi : m_R m_Ksi = m_Ksi_thilde \n",
"\n",
"On obtient m_Ksi_tilde 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 m_Phi, ce qui nous donne une nouvelle matrice des OMs m_Phi_tilde. Pour cela on créer des matrices intérmédiaires:\n",
" - une matrice ( m_Psi ) n par n avec tous ses éléments nuls sauf les colonnes qui contiennent les OMs i et j \n",
" - une matrice ( m_Psi_tilde ) n par n avec tous ses éléments nuls sauf les colonnes qui contiennent les OMs i~ et j~\n",
" - une matrice ( m_interm ) n par n où l'on a soustrait m_Psi à m_Phi pour créer des 0 sur les colonnes des OMs i et j \n",
"\n",
"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 m_Phi_tilde.\n",
"\n",
2020-05-06 15:39:23 +02:00
"On repart de cette nouvelle matrice m_Phi_tilde 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.\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])$"
2020-04-28 16:15:23 +02:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculs : Localisation de Edminston"
]
},
2020-05-07 18:00:34 +02:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fonctions"
2020-04-28 16:15:23 +02:00
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-07 18:00:34 +02:00
"(* Définitions de base nécessaire pour la suite *)\n",
2020-04-29 17:21:24 +02:00
"let ee_ints = AOBasis.ee_ints ao_basis;;\n",
"let m_C = MOBasis.mo_coef mo_basis;;\n",
"let n_ao = Mat.dim1 m_C ;;\n",
2020-05-04 16:14:06 +02:00
"let n_mo = Mat.dim2 m_C ;;\n",
2020-05-07 23:45:41 +02:00
"let multipoles = \n",
" AOBasis.multipole ao_basis;;\n",
2020-05-14 15:54:44 +02:00
" \n",
2020-04-29 17:21:24 +02:00
"let sum a = \n",
" Array.fold_left (fun accu x -> accu +. x) 0. a\n",
2020-05-07 18:00:34 +02:00
" \n",
"type alphaij = {\n",
" alpha_max : float;\n",
" indice_ii : int;\n",
2020-05-18 10:26:28 +02:00
" indice_jj : int;};;\n",
" \n",
"let pi = 3.14159265358979323846264338;;\n"
2020-05-07 18:00:34 +02:00
]
},
2020-05-15 16:31:00 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n",
2020-05-16 11:09:53 +02:00
"(*\n",
2020-05-15 16:31:00 +02:00
"let f_alpha m_C =\n",
"\n",
" let n_mo = Mat.dim2 m_C in\n",
" let t0 = Sys.time () in\n",
" \n",
" let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
" let m_a12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
" \n",
" (* Tableaux temporaires *)\n",
" let m_pqr =\n",
" Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)\n",
" in\n",
" let m_qr_i = Mat.create (n_ao*n_ao) n_mo in\n",
" let m_ri_j = Mat.create (n_ao*n_mo) n_mo in\n",
" let m_ij_k = Mat.create (n_mo*n_mo) n_mo in\n",
" \n",
" Array.iter (fun s ->\n",
" (* Grosse boucle externe sur s *)\n",
" Array.iter (fun r ->\n",
" Array.iter (fun q ->\n",
" Array.iter (fun p ->\n",
" m_pqr.{p,q,r} <- ERI.get_phys ee_ints p q r s\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" (* Conversion d'un tableau a 3 indices en une matrice nao x nao^2 *)\n",
" let m_p_qr =\n",
" Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]\n",
" |> Bigarray.array2_of_genarray\n",
" in\n",
" \n",
" let m_qr_i =\n",
" (* (qr,i) = <i r|q s> = \\sum_p <p r | q s> C_{pi} *)\n",
" gemm ~transa:`T ~c:m_qr_i m_p_qr m_C\n",
" in\n",
" \n",
" let m_q_ri =\n",
" (* Transformation de la matrice (qr,i) en (q,ri) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo)\n",
" in\n",
" \n",
" let m_ri_j =\n",
" (* (ri,j) = <i r | j s> = \\sum_q <i r | q s> C_{bj} *)\n",
" gemm ~transa:`T ~c:m_ri_j m_q_ri m_C\n",
" in\n",
" \n",
" let m_r_ij =\n",
" (* Transformation de la matrice (ri,j) en (r,ij) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo)\n",
" in\n",
" \n",
" let m_ij_k =\n",
" (* (ij,k) = <i k | j s> = \\sum_r <i r | j s> C_{rk} *)\n",
" gemm ~transa:`T ~c:m_ij_k m_r_ij m_C\n",
" in\n",
" \n",
" let m_ijk =\n",
" (* Transformation de la matrice (ei,j) en (e,ij) *)\n",
" Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |]\n",
" |> Bigarray.array3_of_genarray\n",
" in\n",
" \n",
" Array.iter (fun j ->\n",
" Array.iter (fun i ->\n",
" m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});\n",
" m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j} -.\n",
" 0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_C.{s,i} +.\n",
" (m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_C.{s,j})\n",
" ) (Util.array_range 1 n_mo)\n",
" ) (Util.array_range 1 n_mo)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" let t1 = Sys.time () in\n",
" Printf.printf \"t = %f s\\n%!\" (t1 -. t0);\n",
" Mat.init_cols n_mo n_mo ( fun i j ->\n",
" if i= j then 0.\n",
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2. ))))\n",
" );;\n",
"\n",
"(*********************)\n",
"\n",
2020-05-16 11:09:53 +02:00
"f_alpha m_C;;\n",
"*)"
2020-05-15 16:31:00 +02:00
]
},
2020-05-07 18:00:34 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-15 12:33:33 +02:00
"\n",
"\n",
2020-05-07 18:00:34 +02:00
"(*Fonction général de calcul des intégrales*) \n",
2020-04-29 17:21:24 +02:00
"let integral_general g i j =\n",
"Array.map (fun a ->\n",
" let v = \n",
" Array.map (fun b ->\n",
" let u = \n",
" Array.map (fun e ->\n",
" let t = Array.map (fun f ->\n",
" (g a b e f i j) *. ERI.get_phys ee_ints a e b f\n",
" ) (Util.array_range 1 n_ao)\n",
" in sum t\n",
" ) (Util.array_range 1 n_ao)\n",
" in sum u\n",
" ) (Util.array_range 1 n_ao)\n",
" in sum v\n",
") (Util.array_range 1 n_ao)\n",
2020-05-14 15:54:44 +02:00
"|> sum \n",
2020-05-15 12:33:33 +02:00
"\n",
2020-05-15 16:31:00 +02:00
"\n",
2020-05-15 12:33:33 +02:00
"(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n",
"\n",
"let f_alpha m_C =\n",
" let n_mo = Mat.dim2 m_C\n",
" in\n",
" (* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)\n",
" let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> \n",
" integral_general (fun a b e f i j ->\n",
" ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ) *. m_C.{e,i} *. m_C.{f,j}\n",
" ) i j\n",
" )\n",
"\n",
" in\n",
" (* Fonction de calcul de toutes les intégrales A_12 -> Matrice, dépend de m_C *)\n",
" let m_a12 = Mat.init_cols n_mo n_mo (fun i j ->\n",
" integral_general (fun a b e f i j -> m_C.{a,i} *. m_C.{b,j} *. m_C.{e,i} *. m_C.{f,j} \n",
" -. 0.25 *. (( m_C.{e,i} *. m_C.{f,i} -. m_C.{e,j} *. m_C.{f,j} ) \n",
" *. ( m_C.{a,i} *. m_C.{b,i} -. m_C.{a,j} *. m_C.{b,j} ))\n",
" ) i j\n",
" )\n",
"in\n",
"Mat.init_cols n_mo n_mo ( fun i j -> \n",
" if i= j \n",
" then 0. \n",
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))));;\n",
"\n",
"(*********************)\n",
"\n",
"f_alpha m_C;; \n",
2020-05-15 16:31:00 +02:00
"\n"
2020-05-15 12:33:33 +02:00
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n",
"let f_alpha m_C =\n",
"\n",
" let n_mo = Mat.dim2 m_C in\n",
" let t0 = Sys.time () in\n",
" \n",
" let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
" let m_a12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
2020-05-15 17:14:02 +02:00
" let v_d = Vec.init n_mo (fun i -> 0.) in\n",
2020-05-15 12:33:33 +02:00
" \n",
" (* Tableaux temporaires *)\n",
" let m_pqr =\n",
" Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)\n",
" in\n",
" let m_qr_i = Mat.create (n_ao*n_ao) n_mo in\n",
" let m_ri_j = Mat.create (n_ao*n_mo) n_mo in\n",
" let m_ij_k = Mat.create (n_mo*n_mo) n_mo in\n",
" \n",
" Array.iter (fun s ->\n",
" (* Grosse boucle externe sur s *)\n",
" Array.iter (fun r ->\n",
" Array.iter (fun q ->\n",
" Array.iter (fun p ->\n",
" m_pqr.{p,q,r} <- ERI.get_phys ee_ints p q r s\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" (* Conversion d'un tableau a 3 indices en une matrice nao x nao^2 *)\n",
" let m_p_qr =\n",
" Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]\n",
" |> Bigarray.array2_of_genarray\n",
" in\n",
" \n",
" let m_qr_i =\n",
" (* (qr,i) = <i r|q s> = \\sum_p <p r | q s> C_{pi} *)\n",
" gemm ~transa:`T ~c:m_qr_i m_p_qr m_C\n",
" in\n",
" \n",
" let m_q_ri =\n",
" (* Transformation de la matrice (qr,i) en (q,ri) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_qr_i) n_ao (n_ao*n_mo)\n",
" in\n",
" \n",
" let m_ri_j =\n",
" (* (ri,j) = <i r | j s> = \\sum_q <i r | q s> C_{bj} *)\n",
" gemm ~transa:`T ~c:m_ri_j m_q_ri m_C\n",
" in\n",
" \n",
" let m_r_ij =\n",
" (* Transformation de la matrice (ri,j) en (r,ij) *)\n",
" Bigarray.reshape_2 (Bigarray.genarray_of_array2 m_ri_j) n_ao (n_mo*n_mo)\n",
" in\n",
" \n",
" let m_ij_k =\n",
" (* (ij,k) = <i k | j s> = \\sum_r <i r | j s> C_{rk} *)\n",
" gemm ~transa:`T ~c:m_ij_k m_r_ij m_C\n",
" in\n",
" \n",
" let m_ijk =\n",
" (* Transformation de la matrice (ei,j) en (e,ij) *)\n",
" Bigarray.reshape (Bigarray.genarray_of_array2 m_ij_k) [| n_mo ; n_mo ; n_mo |]\n",
" |> Bigarray.array3_of_genarray\n",
" in\n",
" \n",
" Array.iter (fun j ->\n",
" Array.iter (fun i ->\n",
" m_b12.{i,j} <- m_b12.{i,j} +. m_C.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});\n",
" m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_C.{s,j} -.\n",
" 0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_C.{s,i} +.\n",
2020-05-15 17:14:02 +02:00
" (m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_C.{s,j})\n",
" ) (Util.array_range 1 n_mo);\n",
" v_d.{j} <- v_d.{j} +. m_ijk.{j,j,j} *. m_C.{s,j}\n",
2020-05-15 12:33:33 +02:00
" ) (Util.array_range 1 n_mo)\n",
" ) (Util.array_range 1 n_ao);\n",
" \n",
" let t1 = Sys.time () in\n",
" Printf.printf \"t = %f s\\n%!\" (t1 -. t0);\n",
2020-05-15 16:31:00 +02:00
" (Mat.init_cols n_mo n_mo ( fun i j ->\n",
2020-05-15 12:33:33 +02:00
" if i= j then 0.\n",
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2. ))))\n",
2020-05-16 11:09:53 +02:00
" ),Vec.sum v_d);;\n",
2020-05-15 12:33:33 +02:00
"\n",
"(*********************)\n",
"\n",
2020-05-15 16:31:00 +02:00
"f_alpha m_C;;\n",
2020-05-16 11:09:53 +02:00
"(*\n",
"let m_alpha , s_D = f_alpha m_C;;\n",
2020-05-15 16:31:00 +02:00
"\n",
2020-05-16 11:09:53 +02:00
"*)\n",
2020-05-15 16:31:00 +02:00
"\n",
"\n",
2020-05-16 11:09:53 +02:00
"\n"
2020-05-14 15:54:44 +02:00
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-15 12:33:33 +02:00
"(*\n",
2020-05-14 15:54:44 +02:00
"(* Fonction de calcul de tous les alpha ER -> Matrice, dépend de m_a12, m_b12 qui dépendent de m_C *)\n",
2020-05-14 16:06:31 +02:00
"let f_alpha m_C eps=\n",
2020-05-14 15:54:44 +02:00
" let n_mo = Mat.dim2 m_C\n",
" in\n",
" (* Fonction de calcul de toutes les intégrales B_12 -> Matrice, dépend de m_C *)\n",
" let m_b12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
" let m_a12 = Mat.init_cols n_mo n_mo (fun i j -> 0.) in\n",
" Array.iter (fun a ->\n",
" Array.iter (fun b ->\n",
" let mca = Vec.init n_mo (fun i -> m_C.{a,i} *. m_C.{b,i}) in\n",
" Array.iter (fun e ->\n",
" Array.iter (fun f ->\n",
2020-05-14 16:06:31 +02:00
" let integral = ERI.get_phys ee_ints f b e a in\n",
" if abs_float integral > eps then\n",
2020-05-14 15:54:44 +02:00
" Array.iter ( fun i -> \n",
" let mcei = m_C.{e,i} *. m_C.{f,i} in\n",
" Array.iter ( fun j -> \n",
" let x = m_C.{e,i} *. m_C.{f,j} *. integral in\n",
" let mcaij = ( mca.{i} -. mca.{j} ) in\n",
" m_b12.{i,j} <- m_b12.{i,j} +. mcaij *. x;\n",
" m_a12.{i,j} <- m_a12.{i,j} +. m_C.{a,i} *. m_C.{b,j} *. x\n",
" -. 0.25 *. ( mcei -. m_C.{e,j} *. m_C.{f,j} ) *. mcaij *. integral \n",
" ) (Util.array_range 1 n_mo)\n",
" ) (Util.array_range 1 n_mo)\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao)\n",
" ) (Util.array_range 1 n_ao);\n",
"\n",
" Mat.init_cols n_mo n_mo ( fun i j -> \n",
" if i= j then 0. \n",
" else 0.25 *. (acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.)))));;\n",
"\n",
"(*********************)\n",
"\n",
2020-05-15 12:33:33 +02:00
"f_alpha m_C 1.e-5;; \n",
"*)"
2020-05-07 18:00:34 +02:00
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-16 11:09:53 +02:00
"(*\n",
2020-05-07 23:45:41 +02:00
"(* Calcul de D -> critère à maximiser dans ER*)\n",
2020-05-07 18:00:34 +02:00
"let s_D m_C = \n",
" let v_D = \n",
2020-05-14 15:08:25 +02:00
" let n_mo = Mat.dim2 m_C\n",
" in\n",
2020-05-13 15:46:48 +02:00
" let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
2020-05-07 18:00:34 +02:00
" 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",
2020-05-05 15:52:02 +02:00
" ) i j\n",
" )\n",
2020-05-13 15:46:48 +02:00
" in Vec.init n_mo ( fun i -> m_D.{i,i} )\n",
2020-05-07 18:00:34 +02:00
" in Vec.sum v_D ;;\n",
2020-04-29 17:21:24 +02:00
"\n",
2020-05-08 10:49:21 +02:00
"(******************)\n",
2020-05-15 16:31:00 +02:00
"let m_D = Mat.init_cols n_mo n_mo (fun i j ->\n",
" 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",
"let toto = s_D m_C;;\n",
2020-05-16 11:20:50 +02:00
"\n",
2020-05-16 11:09:53 +02:00
"*)"
2020-05-07 23:45:41 +02:00
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Fonction calcul alpha Boys *)\n",
2020-05-14 15:08:25 +02:00
"let f_alpha_boys m_C = \n",
" let n_mo = Mat.dim2 m_C\n",
" in\n",
2020-05-07 23:45:41 +02:00
" let phi_x_phi =\n",
" Multipole.matrix_x multipoles \n",
2020-05-14 15:08:25 +02:00
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
2020-05-07 23:45:41 +02:00
" in \n",
" let phi_y_phi =\n",
" Multipole.matrix_y multipoles \n",
2020-05-14 15:08:25 +02:00
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
2020-05-07 23:45:41 +02:00
" in\n",
" let phi_z_phi =\n",
2020-05-08 10:49:21 +02:00
" Multipole.matrix_z multipoles \n",
2020-05-14 15:08:25 +02:00
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
2020-05-07 23:45:41 +02:00
" in \n",
"\n",
" let m_b12= \n",
2020-05-13 15:46:48 +02:00
" let b12 g = Mat.init_cols n_mo n_mo ( fun i j ->\n",
2020-05-07 23:45:41 +02:00
" (g.{i,i} -. g.{j,j}) *. g.{i,j})\n",
"\n",
" in \n",
" Mat.add (b12 phi_x_phi) ( Mat.add (b12 phi_y_phi) (b12 phi_z_phi))\n",
" in \n",
" let m_a12 =\n",
2020-05-13 15:46:48 +02:00
" let a12 g = Mat.init_cols n_mo n_mo ( fun i j -> \n",
2020-05-11 18:22:34 +02:00
" g.{i,j} *. g.{i,j} -. 0.25 *. ((g.{i,i} -. g.{j,j}) *. (g.{i,i} -. g.{j,j})))\n",
2020-05-07 23:45:41 +02:00
" in\n",
" Mat.add (a12 phi_x_phi) ( Mat.add (a12 phi_y_phi) (a12 phi_z_phi))\n",
" in\n",
2020-05-18 10:26:28 +02:00
" (Mat.init_cols n_mo n_mo ( fun i j -> \n",
2020-05-11 15:25:39 +02:00
" if i=j \n",
" then 0.\n",
2020-05-18 10:26:28 +02:00
" else 0.25 *. acos(-. m_a12.{i,j} /. sqrt((m_a12.{i,j}**2.) +. (m_b12.{i,j}**2.) ))),\n",
" Vec.sum(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",
2020-05-07 23:45:41 +02:00
"\n",
2020-05-08 10:49:21 +02:00
"(****************************)\n",
2020-05-18 10:26:28 +02:00
"\n",
"f_alpha_boys m_C;;\n"
2020-05-07 23:45:41 +02:00
]
},
2020-05-08 10:49:21 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Fonction de calcul de D Boys *)\n",
2020-05-18 10:26:28 +02:00
"(*\n",
2020-05-14 15:08:25 +02:00
"let d_boys m_C = \n",
2020-05-08 10:49:21 +02:00
"\n",
2020-05-14 15:08:25 +02:00
" let phi_x_phi =\n",
" Multipole.matrix_x multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C \n",
2020-05-08 10:49:21 +02:00
" in\n",
" \n",
2020-05-14 15:08:25 +02:00
" (*Util.debug_matrix \"phi_x_phi\" phi_x_phi;*)\n",
2020-05-08 10:49:21 +02:00
" \n",
2020-05-14 15:08:25 +02:00
" let phi_y_phi =\n",
" Multipole.matrix_y multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C \n",
2020-05-08 10:49:21 +02:00
" in\n",
" \n",
2020-05-14 15:08:25 +02:00
" (*Util.debug_matrix \"phi_y_phi\" phi_y_phi;*)\n",
2020-05-08 10:49:21 +02:00
" \n",
2020-05-14 15:08:25 +02:00
" let phi_z_phi =\n",
" Multipole.matrix_z multipoles \n",
" |> MOBasis.mo_matrix_of_ao_matrix ~mo_coef:m_C\n",
2020-05-08 10:49:21 +02:00
"\n",
" in\n",
" \n",
2020-05-14 15:08:25 +02:00
" (*Util.debug_matrix \"phi_z_phi\" phi_z_phi;*)\n",
2020-05-08 10:49:21 +02:00
" \n",
" let v_D_boys = \n",
2020-05-14 15:08:25 +02:00
" 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",
2020-05-08 10:49:21 +02:00
"in\n",
"Vec.sum v_D_boys;;\n",
2020-05-18 10:26:28 +02:00
"*)\n",
2020-05-08 10:49:21 +02:00
"\n",
2020-05-14 08:54:39 +02:00
"(*************************)\n",
2020-05-18 10:26:28 +02:00
"(*\n",
"Multipole.matrix_x multipoles;;\n",
"d_boys m_C;;\n",
"*)"
2020-05-08 10:49:21 +02:00
]
},
2020-05-11 10:58:04 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Test méthode de calcul de alpha et de D ensemble *)\n",
"type alphad = {\n",
"m_alpha : Mat.t;\n",
"d : float;\n",
"}\n",
"\n",
2020-05-11 15:25:39 +02:00
"let m_alpha_d methode m_C = \n",
2020-05-11 10:58:04 +02:00
"\n",
2020-05-18 10:26:28 +02:00
" let alpha_boys , d_boys = f_alpha_boys m_C\n",
2020-05-13 12:34:27 +02:00
" in\n",
2020-05-18 10:26:28 +02:00
" (*let d_boys = d_boys m_C\n",
" in*)\n",
2020-05-16 11:09:53 +02:00
" let alpha_er , d_er = f_alpha m_C\n",
2020-05-13 12:34:27 +02:00
" in\n",
2020-05-16 11:09:53 +02:00
" (*let alpha_er = mat_alpha m_C\n",
2020-05-13 12:34:27 +02:00
" in\n",
2020-05-16 11:09:53 +02:00
" let d_er = s_D m_C\n",
" in*)\n",
2020-05-13 12:34:27 +02:00
" let alpha methode =\n",
" match methode with \n",
" | \"Boys\"\n",
" | \"boys\" -> {m_alpha = alpha_boys; d = d_boys}\n",
" | \"ER\"\n",
" | \"er\" -> {m_alpha = alpha_er; d = d_er}\n",
" | _ -> invalid_arg \"Unknown method, please enter Boys or ER\"\n",
2020-05-11 10:58:04 +02:00
"\n",
"in \n",
"alpha methode;;\n",
"\n",
"(*************************)\n",
"(*\n",
"m_alpha_d \"ER\" ;;\n",
"m_alpha_d \"Boys\" ;;\n",
"\n",
"let methode = \"ER\";;\n",
"\n",
2020-05-14 08:54:39 +02:00
"let alphad = m_alpha_d methode m_C;;\n",
2020-05-11 10:58:04 +02:00
"\n",
"let m_alpha = alphad.m_alpha;;\n",
"let d = alphad.d;;\n",
"*)"
2020-05-07 18:00:34 +02:00
]
},
2020-05-13 12:34:27 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Test norme de alpha *)\n",
"\n",
"let norme m = \n",
" let vec_m = Mat.as_vec m\n",
" in\n",
" let vec2 = Vec.sqr vec_m\n",
"in sqrt(Vec.sum vec2);;\n",
" \n",
"(*\n",
"norme_alpha m_alpha;;\n",
"*)"
]
},
2020-05-07 18:00:34 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Détermination alpha_max et ses indices i et j.\n",
"Si alpha max > pi/2 on soustrait pi/2 à la matrice des alphas de manière récursive *)\n",
2020-05-15 12:33:33 +02:00
"let rec new_m_alpha m_alpha m_C n_rec_alpha=\n",
"\n",
" let n_mo = Mat.dim2 m_C\n",
2020-05-14 15:08:25 +02:00
" in\n",
2020-05-07 18:00:34 +02:00
" let alpha_m =\n",
2020-05-05 15:52:02 +02:00
" \n",
2020-05-07 18:00:34 +02:00
" (*Printf.printf \"%i\\n%!\" n_rec_alpha;*)\n",
2020-05-06 15:39:23 +02:00
" \n",
2020-05-07 18:00:34 +02:00
" if n_rec_alpha == 0 \n",
" then m_alpha \n",
2020-05-13 15:46:48 +02:00
" else Mat.init_cols n_mo n_mo (fun i j -> \n",
2020-05-18 10:26:28 +02:00
" if (m_alpha.{i,j}) > (pi /. 2.) \n",
" then (m_alpha.{i,j} -. ( pi /. 2.))\n",
" else if m_alpha.{i,j} < -. pi /. 2.\n",
" then (m_alpha.{i,j} +. ( pi /. 2.))\n",
2020-05-07 18:00:34 +02:00
" else if m_alpha.{i,j} < 0. \n",
" then -. m_alpha.{i,j}\n",
" else m_alpha.{i,j} )\n",
" in \n",
2020-05-06 15:39:23 +02:00
" \n",
2020-05-07 18:00:34 +02:00
" (*Util.debug_matrix \"alpha_m\" alpha_m;*)\n",
2020-05-05 15:52:02 +02:00
" \n",
2020-05-07 18:00:34 +02:00
" (* Détermination de l'emplacement du alpha max *)\n",
" let max_element3 alpha_m = \n",
" Mat.as_vec alpha_m\n",
" |> iamax\n",
" in\n",
2020-05-13 23:10:33 +02:00
" \n",
2020-05-07 18:00:34 +02:00
" (* indice i du alpha max *)\n",
" let indice_ii = \n",
" let max = max_element3 alpha_m (* Fonction -> constante *)\n",
" in\n",
2020-05-06 15:39:23 +02:00
" \n",
2020-05-07 18:00:34 +02:00
" (*Printf.printf \"%i\\n%!\" max;*)\n",
" \n",
2020-05-13 23:10:33 +02:00
" (max - 1) mod n_mo +1 \n",
2020-05-06 15:39:23 +02:00
" in\n",
2020-05-07 18:00:34 +02:00
" \n",
" (* indice j du alpha max *)\n",
" let indice_jj = \n",
" let max = max_element3 alpha_m (* Fonction -> constante *)\n",
" in\n",
2020-05-13 15:46:48 +02:00
" (max - 1) / n_mo +1\n",
2020-05-06 15:39:23 +02:00
" in\n",
2020-05-07 18:00:34 +02:00
" \n",
" (* Valeur du alpha max*)\n",
" let alpha alpha_m = \n",
" let i = indice_ii \n",
2020-05-04 15:22:13 +02:00
" in\n",
2020-05-07 18:00:34 +02:00
" let j = indice_jj \n",
2020-05-06 15:39:23 +02:00
" in\n",
" \n",
2020-05-07 18:00:34 +02:00
" (*Printf.printf \"%i %i\\n%!\" i j;*)\n",
" \n",
" alpha_m.{i,j}\n",
" \n",
" in\n",
" let alpha_max = alpha alpha_m (* Fonction -> constante *)\n",
2020-04-29 17:21:24 +02:00
" in\n",
2020-05-05 15:52:02 +02:00
" \n",
2020-05-07 18:00:34 +02:00
" (*Printf.printf \"%f\\n%!\" alpha_max;*)\n",
2020-05-05 15:52:02 +02:00
" \n",
2020-05-18 10:26:28 +02:00
" if alpha_max < pi /. 2.\n",
2020-05-07 18:00:34 +02:00
" then {alpha_max; indice_ii; indice_jj}\n",
2020-05-15 12:33:33 +02:00
" else new_m_alpha alpha_m m_C (n_rec_alpha-1);;\n",
2020-05-07 18:00:34 +02:00
"\n",
2020-05-08 10:49:21 +02:00
"(*************************)\n",
2020-05-18 10:26:28 +02:00
"\n",
"let m_alpha,d = f_alpha m_C\n",
2020-05-15 16:31:00 +02:00
"let alphaij = new_m_alpha m_alpha m_C 3;;\n",
2020-05-18 10:26:28 +02:00
"alphaij.alpha_max;;\n"
2020-05-07 18:00:34 +02:00
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-18 10:26:28 +02:00
"let alpha_v loc_deloc alphaij = \n",
" let alpha_loc = alphaij.alpha_max\n",
" in\n",
" let alpha_deloc = alphaij.alpha_max +. (pi /. 4.)\n",
" in\n",
" let choice loc_deloc = \n",
" match loc_deloc with\n",
" |\"loc\" -> alpha_loc\n",
" |\"deloc\" -> alpha_deloc\n",
" | _ -> invalid_arg \"Unknown method, please enter loc or deloc\" \n",
"in choice loc_deloc ;;\n",
"\n",
"\n",
2020-05-07 18:00:34 +02:00
"(* Matrice de rotation 2 par 2 *)\n",
2020-05-16 11:09:53 +02:00
"let f_R alpha =\n",
2020-05-15 16:31:00 +02:00
" Mat.init_cols 2 2 (fun i j -> \n",
" if i=j \n",
" then cos alpha\n",
" else if i>j \n",
" then sin alpha \n",
" else -. sin alpha )\n",
2020-05-18 10:26:28 +02:00
";;\n",
2020-05-08 10:49:21 +02:00
"(*************************)\n",
"\n",
2020-05-16 11:09:53 +02:00
"(*\n",
2020-05-07 18:00:34 +02:00
"let alpha = alphaij.alpha_max;; (* Fonction -> constante *) \n",
2020-05-16 11:09:53 +02:00
"f_R alpha;;\n",
2020-05-18 10:26:28 +02:00
"*)\n",
"alpha_v \"deloc\" alphaij;;\n",
"let alpha = (alpha_v \"loc\" alphaij) ;;\n",
"f_R alpha ;;\n"
2020-05-07 18:00:34 +02:00
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(*Uniquement pour pouvoir tester les fonctions après cette cellules*)\n",
"(*\n",
"(* Indice i et j du alpha max après calcul *)\n",
"let indice_i = alphaij.indice_ii;; (* Fonction -> constante *)\n",
" \n",
"let indice_j = alphaij.indice_jj;; (* Fonction -> constante *)\n",
2020-05-05 15:52:02 +02:00
" \n",
" \n",
2020-05-07 18:00:34 +02:00
"let m_R = f_R alpha;; (* Fonction -> constante *) \n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 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)\n",
"pour appliquer R afin d'effectuer la rotation des orbitales *) (* {1,2} -> 1ere ligne, 2e colonne *)\n",
2020-05-15 12:33:33 +02:00
"let f_Ksi indice_i indice_j m_C =\n",
" let n_ao = Mat.dim1 m_C\n",
"in\n",
"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",
2020-04-29 17:21:24 +02:00
"\n",
2020-05-08 10:49:21 +02:00
"(*************************)\n",
"\n",
2020-05-07 18:00:34 +02:00
"(*\n",
"let m_Ksi = f_Ksi indice_i indice_j m_C;; (* Fonction -> constante *)\n",
2020-05-15 12:33:33 +02:00
"*)\n"
2020-05-07 18:00:34 +02:00
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 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",
2020-05-15 12:33:33 +02:00
"let f_Ksi_tilde m_R m_Ksi m_C = gemm m_Ksi m_R;;\n",
2020-05-08 10:49:21 +02:00
"\n",
"(*************************)\n",
"\n",
2020-05-07 18:00:34 +02:00
"(*\n",
"let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi;; (* Fonction -> constante *)\n",
2020-05-15 16:31:00 +02:00
"*)\n"
2020-05-07 18:00:34 +02:00
]
},
{
"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",
2020-05-06 22:34:25 +02:00
" \n",
2020-05-07 18:00:34 +02:00
"(* Matrice intérmédiare pour l'injection de ksi~ (i~ et j~) dans la matrice Phi *)\n",
2020-05-18 10:26:28 +02:00
"(*\n",
2020-05-15 12:33:33 +02:00
"let f_Psi_tilde m_Ksi_tilde indice_i indice_j m_C=\n",
" let n_mo = Mat.dim2 m_C\n",
" in\n",
" let n_ao = Mat.dim1 m_C\n",
2020-05-18 10:26:28 +02:00
"in\n",
2020-05-15 12:33:33 +02:00
"Mat.init_cols n_ao n_mo (fun i j -> \n",
2020-05-07 18:00:34 +02:00
" if j=indice_i \n",
" then m_Ksi_tilde.{i,1}\n",
2020-05-18 10:26:28 +02:00
" else if j=indice_j \n",
" then m_Ksi_tilde.{i,2}\n",
2020-05-07 18:00:34 +02:00
" else 0.)\n",
2020-05-18 10:26:28 +02:00
" \n",
2020-05-07 18:00:34 +02:00
";;\n",
2020-05-18 10:26:28 +02:00
"*)\n",
2020-05-07 18:00:34 +02:00
"(* Matrice intermédiaire pour supprimer ksi (i et j) dans la matrice Phi *) \n",
2020-05-18 10:26:28 +02:00
"(*\n",
2020-05-15 12:33:33 +02:00
"let f_Psi m_Ksi indice_i indice_j m_C = \n",
" let n_mo = Mat.dim2 m_C\n",
" in\n",
" let n_ao = Mat.dim1 m_C\n",
2020-05-18 10:26:28 +02:00
"in\n",
2020-05-15 12:33:33 +02:00
"Mat.init_cols n_ao n_mo (fun i j -> \n",
2020-05-07 18:00:34 +02:00
" if j=indice_i \n",
" then m_Ksi.{i,1}\n",
" else if j=indice_j \n",
" then m_Ksi.{i,2}\n",
" else 0.)\n",
";;\n",
2020-05-18 10:26:28 +02:00
"*)\n",
"let f_k mat indice_i indice_j m_C = \n",
"let n_mo = Mat.dim2 m_C\n",
" in\n",
" let n_ao = Mat.dim1 m_C\n",
"in\n",
"Mat.init_cols n_ao n_mo (fun i j -> \n",
" if j=indice_i \n",
" then mat.{i,1}\n",
" else if j=indice_j \n",
" then mat.{i,2}\n",
" else 0.)\n",
2020-05-08 10:49:21 +02:00
"(*************************)\n",
"\n",
2020-05-07 18:00:34 +02:00
"(*\n",
"let m_Psi = f_Psi m_Ksi indice_i indice_j;; (* Fonction -> constante *)\n",
"\n",
"let m_Psi_tilde = f_Psi_tilde m_Ksi_tilde indice_i indice_j;; (* Fonction -> constante *)\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* 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",
2020-05-08 10:49:21 +02:00
"\n",
"(*************************)\n",
"\n",
2020-05-06 22:34:25 +02:00
"(*\n",
2020-05-07 18:00:34 +02:00
"let m_interm = f_interm m_C m_Psi;; (* Fonction -> constante *)\n",
"\n",
"let new_m_C m_C= Mat.add m_Psi_tilde m_interm;;\n",
"\n",
"let m_new_m_C = new_m_C m_C;;\n",
2020-05-06 22:34:25 +02:00
"*)"
]
},
2020-05-14 08:54:39 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-15 16:31:00 +02:00
"(* Test localisation matrice rectangulaire et partielle *)\n",
"(*let toto = [4];;\n",
2020-05-15 12:33:33 +02:00
"let occ_m_C m_C toto= Mat.init_cols 4 3 (fun i j ->\n",
2020-05-14 08:54:39 +02:00
" if not (List.mem j toto) \n",
" then m_C.{i,j}\n",
" else 0.);;\n",
" \n",
2020-05-15 16:31:00 +02:00
"let occ = occ_m_C m_C toto;;\n",
"*)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-18 10:26:28 +02:00
"(*Fonction de création d'une list d'entier à partir d'un vecteur de float*)\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",
2020-05-15 16:31:00 +02:00
"(* Fonction créant une liste à partir des éléments manquant d'une autre liste, dans l'intervalle [1 ; n_mo] *)\n",
"let miss_elem mat list = \n",
" let n_mo = Mat.dim2 mat\n",
" in\n",
" let vec = Vec.init (n_mo) (fun i ->\n",
" if List.mem i list\n",
" then 0.\n",
" else float_of_int(i))\n",
" in\n",
2020-05-18 10:26:28 +02:00
" let list_int = int_list vec \n",
"in\n",
"List.filter (fun x -> x > 0) list_int;;\n",
2020-05-15 16:31:00 +02:00
"\n",
"(* Fonction de séparation d'une matrice en 2 sous matrice, la première matrice correspondant aux colonnes de la matrice dont le numéro est présent\n",
"dans la liste et la seconde à celles dont le numéro de colonne n'est pas présent dans la liste *)\n",
"let split_mat mat list =\n",
" let vec_of_mat = Mat.to_col_vecs mat\n",
" in\n",
" let f a = vec_of_mat.(a-1)\n",
" in\n",
" let vec_list_1 = List.map f list\n",
" in\n",
" let list_2 = miss_elem mat list\n",
" in\n",
" let vec_list_2 = List.map f list_2\n",
"in (Mat.of_col_vecs_list vec_list_1,Mat.of_col_vecs_list vec_list_2);;\n",
"\n",
"m_C;;\n",
"\n",
2020-05-15 17:14:02 +02:00
"let list_om = [1;3]\n",
2020-05-15 16:31:00 +02:00
"\n",
2020-05-18 10:26:28 +02:00
"let m_occ , m_vir = split_mat m_C list_om;;\n",
"\n",
"(* Liste des OMs occupées *)\n",
"let list_occ = \n",
" let vec_occ = Vec.init (nocc) (fun i -> float_of_int(i))\n",
"in int_list vec_occ;;\n"
2020-05-14 08:54:39 +02:00
]
},
2020-05-11 10:58:04 +02:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Boucle localisation avec choix de la méthode de localisation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(* Localisation de Edminstion ou de Boys *)\n",
2020-05-13 12:34:27 +02:00
"\n",
2020-05-11 10:58:04 +02:00
"(* Calcul de la nouvelle matrice des coefficient après n rotation d'orbitales *)\n",
2020-05-18 10:26:28 +02:00
"let rec final_m_C m_C methode loc_deloc epsilon n prev_critere_D cc=\n",
2020-05-11 10:58:04 +02:00
"\n",
2020-05-13 15:46:48 +02:00
" Printf.printf \"%i\\n%!\" n;\n",
2020-05-11 10:58:04 +02:00
"\n",
2020-05-12 11:17:15 +02:00
" (*Util.debug_matrix \"m_C\" m_C;*)\n",
2020-05-11 10:58:04 +02:00
"\n",
"if n == 0 \n",
" then m_C\n",
" else\n",
" \n",
" (* Fonction de calcul de la nouvelle matrice de coef après rotation d'un angle alpha *)\n",
2020-05-18 10:26:28 +02:00
" let new_m_C m_C methode loc_deloc =\n",
2020-05-11 10:58:04 +02:00
" \n",
2020-05-12 10:55:33 +02:00
" (* Fonction de pattern matching en fonction de la méthode *)\n",
2020-05-11 15:25:39 +02:00
" let alphad = m_alpha_d methode m_C \n",
2020-05-11 10:58:04 +02:00
" in\n",
" \n",
" (* D critère à maximiser *)\n",
" let critere_D = alphad.d \n",
" in\n",
2020-05-13 12:34:27 +02:00
" \n",
2020-05-11 10:58:04 +02:00
" Printf.printf \"%f\\n%!\" critere_D;\n",
" \n",
" (* Matrice des alphas *)\n",
" let m_alpha = alphad.m_alpha\n",
" in\n",
2020-05-13 12:34:27 +02:00
" let norme_alpha = norme m_alpha\n",
" in\n",
" \n",
2020-05-13 15:46:48 +02:00
" Printf.printf \"%f\\n%!\" norme_alpha;\n",
2020-05-13 12:34:27 +02:00
" \n",
2020-05-12 11:17:15 +02:00
" (*Util.debug_matrix \"m_alpha\" m_alpha;*)\n",
2020-05-11 10:58:04 +02:00
"\n",
" (* alphaij contient le alpha max ainsi que ses indices i et j *)\n",
2020-05-13 12:34:27 +02:00
" let n_rec_alpha = 10 (* Nombre ditération max pour réduire les valeurs de alpha *)\n",
" in\n",
2020-05-15 12:33:33 +02:00
" let alphaij = new_m_alpha m_alpha m_C n_rec_alpha (* Fonction -> constante *)\n",
2020-05-11 10:58:04 +02:00
" in\n",
"\n",
2020-05-12 11:17:15 +02:00
" (* Valeur de alpha max après calcul *) (* Epsilon = Pas <1. , 1. -> normal, sinon Pas plus petit *)\n",
2020-05-18 10:26:28 +02:00
" let alpha = (alpha_v loc_deloc alphaij) *. epsilon (* Fonction -> constante *)\n",
2020-05-11 10:58:04 +02:00
" in\n",
"\n",
2020-05-18 10:26:28 +02:00
" Printf.printf \"%f\\n%!\" alpha;\n",
2020-05-15 16:31:00 +02:00
" \n",
2020-05-11 10:58:04 +02:00
" (* Indice i et j du alpha max après calcul *)\n",
" let indice_i = alphaij.indice_ii (* Fonction -> constante *)\n",
" in\n",
" let indice_j = alphaij.indice_jj (* Fonction -> constante *)\n",
" in\n",
"\n",
2020-05-13 12:34:27 +02:00
" (*Printf.printf \"%i %i\\n%!\" indice_i indice_j;*)\n",
2020-05-15 16:31:00 +02:00
" \n",
" (* Matrice de rotation *)\n",
2020-05-16 11:09:53 +02:00
" let m_R = f_R alpha (* Fonction -> constante *)\n",
2020-05-11 10:58:04 +02:00
" in\n",
"\n",
2020-05-16 11:09:53 +02:00
" (*Util.debug_matrix \"m_R\" m_R;*)\n",
2020-05-11 10:58:04 +02:00
"\n",
" (* Matrice qui va subir la rotation *)\n",
" let m_Ksi = f_Ksi indice_i indice_j m_C (* Fonction -> constante *)\n",
" in\n",
"\n",
2020-05-12 11:17:15 +02:00
" (*Util.debug_matrix \"m_Ksi\" m_Ksi;*)\n",
2020-05-11 10:58:04 +02:00
"\n",
" (* Matrice ayant subit la rotation *)\n",
2020-05-15 12:33:33 +02:00
" let m_Ksi_tilde = f_Ksi_tilde m_R m_Ksi m_C (* Fonction -> constante *)\n",
2020-05-11 10:58:04 +02:00
" in\n",
"\n",
2020-05-12 11:17:15 +02:00
" (*Util.debug_matrix \"m_Ksi_tilde\" m_Ksi_tilde;*)\n",
2020-05-11 10:58:04 +02:00
"\n",
" (* Matrice pour supprimerles coef des orbitales i et j dans la matrice des coef *)\n",
2020-05-18 10:26:28 +02:00
" let m_Psi = f_k m_Ksi indice_i indice_j m_C(* Fonction -> constante *)\n",
2020-05-11 10:58:04 +02:00
" in\n",
"\n",
2020-05-12 11:17:15 +02:00
" (*Util.debug_matrix \"m_Psi\" m_Psi;*)\n",
2020-05-11 10:58:04 +02:00
"\n",
" (* Matrice pour ajouter les coef des orbitales i~ et j~ dans la matrice des coef *)\n",
2020-05-18 10:26:28 +02:00
" let m_Psi_tilde = f_k m_Ksi_tilde indice_i indice_j m_C (* Fonction -> constante *)\n",
2020-05-11 10:58:04 +02:00
" in\n",
"\n",
2020-05-12 11:17:15 +02:00
" (*Util.debug_matrix \"m_Psi_tilde\" m_Psi_tilde;*)\n",
2020-05-11 10:58:04 +02:00
"\n",
" (* Matrice avec les coef des orbitales i et j remplacés par 0 *)\n",
" let m_interm = f_interm m_C m_Psi (* Fonction -> constante *)\n",
" in\n",
2020-05-13 23:10:33 +02:00
" \n",
2020-05-12 11:17:15 +02:00
" (*Util.debug_matrix \"m_interm\" m_interm;*)\n",
2020-05-11 10:58:04 +02:00
" \n",
" (* Matrice après rotation *)\n",
2020-05-13 15:46:48 +02:00
" ( Mat.add m_Psi_tilde m_interm, critere_D)\n",
2020-05-11 10:58:04 +02:00
" in\n",
2020-05-18 10:26:28 +02:00
" let m_new_m_C , critere_D = new_m_C m_C methode loc_deloc (* Fonction -> constante *)\n",
2020-05-13 12:34:27 +02:00
" in\n",
2020-05-14 08:54:39 +02:00
" let diff = prev_critere_D -. critere_D +. 1.\n",
2020-05-13 12:34:27 +02:00
" \n",
2020-05-11 10:58:04 +02:00
"in\n",
2020-05-13 15:46:48 +02:00
"\n",
2020-05-12 11:17:15 +02:00
"(*Util.debug_matrix \"new_alpha_m\" (f_alpha m_C);*)\n",
2020-05-13 12:34:27 +02:00
"(*Util.debug_matrix \"m_new_m_C\" m_new_m_C;*)\n",
2020-05-11 10:58:04 +02:00
"\n",
2020-05-13 15:46:48 +02:00
"if diff**2. < cc**2.\n",
2020-05-13 12:34:27 +02:00
" then m_new_m_C\n",
" else\n",
"\n",
2020-05-18 10:26:28 +02:00
"final_m_C m_new_m_C methode loc_deloc epsilon (n-1) critere_D cc;;"
2020-05-11 15:25:39 +02:00
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-12 11:17:15 +02:00
"(* Calcul *)\n",
2020-05-18 10:26:28 +02:00
"(* Fonction / Matrice des coef / Méthode(\"Boys\" ou \"ER\") / Localisation ou non (\"loc\" ou \"deloc\") / Pas(<=1.) \n",
"/ Nombre d'itérations max / 0. (valeur de D pour initier la boucle) / critère de convergence sur D*)\n",
"let new_m = final_m_C m_C \"boys\" \"loc\" 1. 1 0. 10e-7;;\n"
2020-05-11 10:58:04 +02:00
]
},
2020-05-13 12:34:27 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
2020-05-14 08:54:39 +02:00
"source": [
2020-05-18 10:26:28 +02:00
"final_m_C m_C \"boys\" \"loc\" 1. 1 0. 10e-7;;"
2020-05-05 15:52:02 +02:00
]
},
2020-05-06 22:34:25 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
2020-05-07 18:00:34 +02:00
"source": [
2020-05-18 10:26:28 +02:00
"final_m_C new_m \"Boys\" \"deloc\" 1. 1 0. 10e-7;;\n"
2020-05-07 18:00:34 +02:00
]
2020-05-06 22:34:25 +02:00
},
2020-05-05 15:52:02 +02:00
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
2020-05-18 10:26:28 +02:00
"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",
2020-05-06 22:34:25 +02:00
"\n",
2020-05-18 10:26:28 +02:00
"let int_list vec = \n",
" let float_list = Vec.to_list vec\n",
2020-05-07 14:45:41 +02:00
" in\n",
2020-05-18 10:26:28 +02:00
" let g a = int_of_float(a)\n",
"in List.map g float_list;;\n",
2020-05-07 14:45:41 +02:00
"\n",
2020-05-18 10:26:28 +02:00
"int_list toto;;"
2020-05-07 14:45:41 +02:00
]
},
2020-04-28 16:15:23 +02:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# OLD ->"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"let ints = MOBasis.ee_ints mo_basis in\n",
"ERI.get_phys ints 1 2 1 2;;"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "OCaml 4.10.0",
"language": "OCaml",
"name": "ocaml-jupyter"
},
"language_info": {
"codemirror_mode": "text/x-ocaml",
"file_extension": ".ml",
"mimetype": "text/x-ocaml",
"name": "OCaml",
"nbconverter_exporter": null,
"pygments_lexer": "OCaml",
"version": "4.10.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
