{
 "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
}