From 3c04abd38f5422d75174ef6c1b6db15fb0fc3187 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 19 Apr 2020 23:29:00 +0200 Subject: [PATCH] Localization --- .gitattributes | 2 + Localization.ipynb | 622 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 624 insertions(+) create mode 100644 .gitattributes create mode 100644 Localization.ipynb diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..305a9b0 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,2 @@ +*.ipynb filter=nbstripout +*.ipynb diff=ipynb diff --git a/Localization.ipynb b/Localization.ipynb new file mode 100644 index 0000000..c33f445 --- /dev/null +++ b/Localization.ipynb @@ -0,0 +1,622 @@ +{ + "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": { + "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_{n}" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "heading_collapsed": true + }, + "source": [ + "### Generation du fichier xyz" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hidden": true + }, + "source": [ + "Pour faire plusieurs chaines d'hydrogene, il faut d'abord\n", + "construire une fonction qui genere un fichier xyz pour\n", + "une chaine de `n` hydrogenes avec une distance H-H de `d` bohr.\n", + "Le fichier xyz est en Angstrom, donc il faut utiliser\n", + "`Constants.a0` pour convertir les bohr en Angstroms.\n", + "\n", + "Voila un exemple de ce qu'on doit obtenir pour 20 hydrogenes:\n", + "\n", + "```\n", + " 20\n", + "Hydrogen chain, d=1.8 Bohr\n", + "H -4.286335 0.000000 0.000000\n", + "H -3.333816 0.000000 0.000000\n", + "H -2.381297 0.000000 0.000000\n", + "H -1.428778 0.000000 0.000000\n", + "H -0.476259 0.000000 0.000000\n", + "H 0.476259 0.000000 0.000000\n", + "H 1.428778 0.000000 0.000000\n", + "H 2.381297 0.000000 0.000000\n", + "H 3.333816 0.000000 0.000000\n", + "H 4.286335 0.000000 0.000000\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hidden": true + }, + "outputs": [], + "source": [ + "let string_of_xyz d n =\n", + "\" 10\n", + "Hydrogen chain, d=1.8 Bohr\n", + "H -4.286335 0.000000 0.000000\n", + "H -3.333816 0.000000 0.000000\n", + "H -2.381297 0.000000 0.000000\n", + "H -1.428778 0.000000 0.000000\n", + "H -0.476259 0.000000 0.000000\n", + "H 0.476259 0.000000 0.000000\n", + "H 1.428778 0.000000 0.000000\n", + "H 2.381297 0.000000 0.000000\n", + "H 3.333816 0.000000 0.000000\n", + "H 4.286335 0.000000 0.000000\n", + "\"\n", + "\n", + "let xyz_string =\n", + " string_of_xyz 1.8 10" + ] + }, + { + "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", + "Voila la base STO-6G pour l'hydrogene:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hidden": true + }, + "outputs": [], + "source": [ + "let basis_string = \n", + "\"\n", + "HYDROGEN\n", + "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", + "\n", + "\"" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hidden": true + }, + "source": [ + "Une orbitale atomique STO-6G centree sur l'atome A est composee d'une contraction de 6 Gaussiennes:\n", + "$$\n", + "\\chi(r) = \\sum_{i=1}^6 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 pour generer des orbitales moleculaires" + ] + }, + { + "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) ]\n", + "in\n", + "\n", + "plot data \"test_data.png\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Calcul des integrales necessaires" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "let dx, x_values = \n", + " let n = 101 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) -. 10. ,\n", + " coord.(Array.length coord -1) +. 10.\n", + " in\n", + "\n", + " let dx =\n", + " (xmax -. xmin) /. (float_of_int n -. 1.)\n", + " in\n", + " dx, Array.init n (fun i -> xmin +. (float_of_int i)*.dx)\n", + "in\n", + "\n", + "let dy, y_values = \n", + " let n = 51 in\n", + " \n", + " let ymin, ymax =\n", + " let coord =\n", + " Array.map snd nuclei\n", + " |> Array.map (fun a -> Coordinate.(get Y) a)\n", + " in\n", + " Array.sort compare coord;\n", + " coord.(0) -. 10. ,\n", + " coord.(Array.length coord -1) +. 10.\n", + " in\n", + "\n", + " let dy =\n", + " (ymax -. ymin) /. (float_of_int n -. 1.)\n", + " in\n", + " dy, Array.init n (fun i -> ymin +. (float_of_int i)*.dy)\n", + "in\n", + "\n", + "\n", + "let dz, z_values = \n", + " let n = 51 in\n", + " \n", + " let zmin, zmax =\n", + " let coord =\n", + " Array.map snd nuclei\n", + " |> Array.map (fun a -> Coordinate.(get Z) a)\n", + " in\n", + " Array.sort compare coord;\n", + " coord.(0) -. 10. ,\n", + " coord.(Array.length coord -1) +. 10.\n", + " in\n", + "\n", + " let dz =\n", + " (zmax -. zmin) /. (float_of_int n -. 1.)\n", + " in\n", + " dz, Array.init n (fun i -> zmin +. (float_of_int i)*.dz)\n", + "in\n", + "\n", + "\n", + "let points = \n", + " Array.map (fun x -> \n", + " Array.map (fun y -> \n", + " Array.map (fun z -> \n", + " Coordinate.make_angstrom\n", + " { Coordinate.\n", + " x ; y ; z \n", + " } \n", + " ) x_values\n", + " ) y_values\n", + " ) z_values\n", + " |> Array.map Array.to_list\n", + " |> Array.to_list\n", + " |> List.concat\n", + " |> Array.concat\n", + "in\n", + "\n", + "let dv = dx *. dy *. dz in\n", + "Printf.printf \"%f %f %f %f\\n%!\" dx dy dz dv;\n", + "\n", + "let ao_val = \n", + " Array.map (fun point -> \n", + " AOBasis.values ao_basis point\n", + " ) points\n", + " |> Mat.of_col_vecs\n", + "in\n", + "\n", + "(*\n", + "let r = \n", + " gemm ~transb:`T ao_val ao_val \n", + "in\n", + "Mat.scal dv r;\n", + "r;;\n", + "AOBasis.overlap ao_basis |> Overlap.matrix ;;\n", + "*)\n", + "\n", + "let ao_x_val = \n", + " Array.map (fun point -> \n", + " let r = AOBasis.values ao_basis point in\n", + " scal (Coordinate.(get X) point) r;\n", + " r\n", + " ) points\n", + " |> Mat.of_col_vecs\n", + "in\n", + "\n", + "let r = \n", + " gemm ~transb:`T ao_val ao_x_val \n", + "in\n", + "Mat.scal dv r;\n", + "r;;\n", + "AOBasis.multipole ao_basis |> Multipole.matrix_x ;;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Localisation des orbitales" + ] + }, + { + "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 +}