Hartree-Fock
This commit is contained in:
parent
3c04abd38f
commit
c9b0ef0dfd
778
HartreeFock.ipynb
Normal file
778
HartreeFock.ipynb
Normal file
@ -0,0 +1,778 @@
|
||||
{
|
||||
"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",
|
||||
"#cd \"/home/scemama/QCaml\";;\n",
|
||||
"\n",
|
||||
"#use \"topfind\";;\n",
|
||||
"#require \"jupyter.notebook\";;\n",
|
||||
"\n",
|
||||
"#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": [
|
||||
"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",
|
||||
"\"\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"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": [
|
||||
"let basis_string = \"\n",
|
||||
"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",
|
||||
"\n",
|
||||
"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",
|
||||
"\n",
|
||||
"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
"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",
|
||||
" Simulation.ao_basis simulation"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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 nocc = 5 (* On a 10 electrons, donc 5 orbitales occupees. *)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"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": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"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
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user