{
 "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/ydamour/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\" (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": [
    "let local;;\n"
   ]
  }
 ],
 "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
}