{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test of F12 matrix elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/scemama/qp2/external/opam/4.07.1/lib/bytes: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/base64: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/base64/base64.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/compiler-libs: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/compiler-libs/ocamlcommon.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/result: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/result/result.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/ppx_deriving/runtime: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/ppx_deriving/runtime/ppx_deriving_runtime.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/ppx_deriving_yojson/runtime: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/ppx_deriving_yojson/runtime/ppx_deriving_yojson_runtime.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/unix.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/uuidm: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/uuidm/uuidm.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/easy-format: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/easy-format/easy_format.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/biniou: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/biniou/biniou.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/yojson: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/yojson/yojson.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/jupyter: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/jupyter/jupyter.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/jupyter/notebook: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/jupyter/notebook/jupyter_notebook.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/bigarray.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/lacaml: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/lacaml/lacaml.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/astring: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/astring/astring.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/cmdliner: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/cmdliner/cmdliner.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/seq: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/stdlib-shims: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/stdlib-shims/stdlib_shims.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/fmt: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/fmt/fmt.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/fmt/fmt_cli.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/fmt/fmt_tty.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/alcotest: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/alcotest/alcotest.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/ocaml/str.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/zarith: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/zarith/zarith.cma: loaded\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/getopt: added to search path\n",
      "/home/scemama/qp2/external/opam/4.07.1/lib/getopt/getopt.cma: loaded\n"
     ]
    }
   ],
   "source": [
    "#cd \"/home/scemama/QCaml\";;\n",
    "#use \"topfind\";;\n",
    "#require \"jupyter.notebook\";;\n",
    "\n",
    "#require \"lacaml\";;\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\";;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modules to load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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",
    "open Lacaml.D;;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Printers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "val pp_mo : Format.formatter -> MOBasis.t -> unit = <fun>\n"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#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",
    "let pp_mo ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;\n",
    "#install_printer pp_mo;;\n",
    "#install_printer DeterminantSpace.pp;;\n",
    "#install_printer SpindeterminantSpace.pp;;\n",
    "#install_printer Bitstring.pp;;\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "val basis_filename : string = \"/home/scemama/qp2/data/basis/6-31g\"\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val aux_basis_filename : string = \"/home/scemama/qp2/data/basis/cc-pvdz\"\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val nuclei : Nuclei.t = [|(Element.Be,   0.0000   0.0000   0.0000)|]\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val frozen_core : bool = false\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val multiplicity : int = 1\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val state : int = 1\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val basis : Basis.t =\n",
       "  \n",
       "                          Atomic Basis set\n",
       "                          ----------------\n",
       "\n",
       "-----------------------------------------------------------------------\n",
       " #   Angular   Coordinates (Bohr)       Exponents       Coefficients\n",
       "    Momentum  X        Y        Z\n",
       "-----------------------------------------------------------------------\n",
       "  1-3  S   0.0000   0.0000   0.0000   1.26458570e+03    1.94480000e-03\n",
       "                                      1.89936810e+02    1.48351000e-02\n",
       "                                      4.31590890e+01    7.20906000e-02\n",
       "                                      1.20986630e+01    2.37154200e-01\n",
       "                                      3.80632320e+00    4.69198700e-01\n",
       "                                      1.27289030e+00    3.56520200e-01\n",
       "                                    \n",
       "                                      3.19646310e+00   -1.12648700e-01\n",
       "                                      7.47813300e-01   -2.29506400e-01\n",
       "                                      2.19966300e-01    1.18691670e+00\n",
       "                                    \n",
       "                                      8.23099000e-02    1.00000000e+00\n",
       "                                    \n",
       "                                    \n",
       "-----------------------------------------------------------------------\n",
       "  3-8  P   0.0000   0.0000   0.0000   3.19646310e+00    5.59802000e-02\n",
       "                                      7.47813300e-01    2.61550600e-01\n",
       "                                      2.19966300e-01    7.93972300e-01\n",
       "                                    \n",
       "                                      8.23099000e-02    1.00000000e+00\n",
       "                                    \n",
       "                                    \n",
       "-----------------------------------------------------------------------\n",
       "\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val aux_basis : Basis.t =\n",
       "  \n",
       "                          Atomic Basis set\n",
       "                          ----------------\n",
       "\n",
       "-----------------------------------------------------------------------\n",
       " #   Angular   Coordinates (Bohr)       Exponents       Coefficients\n",
       "    Momentum  X        Y        Z\n",
       "-----------------------------------------------------------------------\n",
       "  1-6  S   0.0000   0.0000   0.0000   1.26458570e+03    1.94480000e-03\n",
       "                                      1.89936810e+02    1.48351000e-02\n",
       "                                      4.31590890e+01    7.20906000e-02\n",
       "                                      1.20986630e+01    2.37154200e-01\n",
       "                                      3.80632320e+00    4.69198700e-01\n",
       "                                      1.27289030e+00    3.56520200e-01\n",
       "                                    \n",
       "                                      3.19646310e+00   -1.12648700e-01\n",
       "                                      7.47813300e-01   -2.29506400e-01\n",
       "                                      2.19966300e-01    1.18691670e+00\n",
       "                                    \n",
       "                                      8.23099000e-02    1.00000000e+00\n",
       "                                    \n",
       "                                      2.94000000e+03    6.80000000e-04\n",
       "                                      4.41200000e+02    5.23600000e-03\n",
       "                                      1.00500000e+02    2.66060000e-02\n",
       "                                      2.84300000e+01    9.99930000e-02\n",
       "                                      9.16900000e+00    2.69702000e-01\n",
       "                                      3.19600000e+00    4.51469000e-01\n",
       "                                      1.15900000e+00    2.95074000e-01\n",
       "                                      1.81100000e-01    1.25870000e-02\n",
       "                                    \n",
       "                                      2.94000000e+03   -1.23000000e-04\n",
       "                                      4.41200000e+02   -9.66000000e-04\n",
       "                                      1.00500000e+02   -4.83100000e-03\n",
       "                                      2.84300000e+01   -1.93140000e-02\n",
       "                                      9.16900000e+00   -5.32800000e-02\n",
       "                                      3.19600000e+00   -1.20723000e-01\n",
       "                                      1.15900000e+00   -1.33435000e-01\n",
       "                                      1.81100000e-01    5.30767000e-01\n",
       "                                    \n",
       "                                      5.89000000e-02    1.00000000e+00\n",
       "                                    \n",
       "                                    \n",
       "-----------------------------------------------------------------------\n",
       "  3-14 P   0.0000   0.0000   0.0000   3.19646310e+00    5.59802000e-02\n",
       "                                      7.47813300e-01    2.61550600e-01\n",
       "                                      2.19966300e-01    7.93972300e-01\n",
       "                                    \n",
       "                                      8.23099000e-02    1.00000000e+00\n",
       "                                    \n",
       "                                      3.61900000e+00    2.91110000e-02\n",
       "                                      7.11000000e-01    1.69365000e-01\n",
       "                                      1.95100000e-01    5.13458000e-01\n",
       "                                    \n",
       "                                      6.01800000e-02    1.00000000e+00\n",
       "                                    \n",
       "                                    \n",
       "-----------------------------------------------------------------------\n",
       " 19-24 D   0.0000   0.0000   0.0000   2.38000000e-01    1.00000000e+00\n",
       "                                    \n",
       "                                    \n",
       "-----------------------------------------------------------------------\n",
       "\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val f12 : F12factor.t =\n",
       "  {F12factor.expo_s = 1.;\n",
       "   gaussian =\n",
       "    {GaussianOperator.coef_g =\n",
       "      [  -0.314400  -0.303700  -0.168100  -0.098110  -0.060240  -0.037260 ];\n",
       "     expo_sg =\n",
       "      [   0.220900   1.004000   3.622000  12.160000  45.870000 254.400000 ];\n",
       "     expo_sg_inv =\n",
       "      [   4.526935   0.996016   0.276091   0.082237   0.021801   0.003931 ]}}\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val charge : int = 0\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val simulation : Simulation.t = <abstr>\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "let basis_filename  = \"/home/scemama/qp2/data/basis/6-31g\" \n",
    "let aux_basis_filename = \"/home/scemama/qp2/data/basis/cc-pvdz\" \n",
    "let nuclei = Nuclei.of_zmt_string \"be\" \n",
    "let frozen_core = false\n",
    "let multiplicity = 1\n",
    "let state = 1\n",
    "\n",
    "let basis =  Basis.of_nuclei_and_basis_filenames ~nuclei [basis_filename] \n",
    "let aux_basis =  Basis.of_nuclei_and_basis_filenames ~nuclei (basis_filename :: aux_basis_filename :: []) \n",
    "let f12 = F12factor.gaussian_geminal 1.0 \n",
    "let charge = 0 \n",
    "\n",
    "\n",
    "let simulation =\n",
    "    Simulation.make \n",
    "        ~f12 ~charge ~multiplicity ~nuclei\n",
    "        ~cartesian:true\n",
    "        basis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hartree-Fock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15 significant shell pairs computed in 0.004974 seconds\n",
      "1\n",
      "2\n",
      "5\n",
      "6\n",
      "Computed ERIs in 0.026928 seconds\n",
      "MOs =\n",
      "\n",
      "\n",
      "          -- 1 --     -- 2 --    -- 3 --   -- 4 --      -- 5 --\n",
      "            1    0.997996    -0.22266          0         0           -0\n",
      "            2   0.0161427    0.288767 2.2404E-15         0 -7.73804E-16\n",
      "            3           0 3.58344E-16   0.258415  0.037609    0.0381979\n",
      "                      ...         ...        ...       ...          ...\n",
      "            7           0           0    0.78727  0.114442     0.116234\n",
      "            8           0           0  -0.128237 0.0798944     0.789475\n",
      "            9           0           0  -0.100874  0.791587   -0.0964853\n",
      "      \n",
      "\n",
      "            -- 6 --      -- 7 --     -- 8 --      -- 9 --\n",
      "            1    0.00338551            0          -0            0\n",
      "            2       2.01724 -3.80748E-15 4.70544E-15 -6.99768E-16\n",
      "            3   1.56069E-15        1.292    0.249105    0.0981672\n",
      "                        ...          ...         ...          ...\n",
      "            7   -2.3622E-15     -1.05652    -0.20374   -0.0802898\n",
      "            8   1.64398E-15     0.172096   -0.528064     -0.92522\n",
      "            9   2.61906E-15     0.135375   -0.918775     0.549573\n",
      "      \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "val hf : HartreeFock.t = \n",
       "======================================================================\n",
       "                        Restricted Hartree-Fock                      \n",
       "======================================================================\n",
       "\n",
       "    ------------------------------------------------------------\n",
       "        #     HF energy       Convergence  HOMO-LUMO\n",
       "    ------------------------------------------------------------\n",
       "        1    -14.35428398      3.4264e-01     0.3789\n",
       "        2    -14.52767889      1.8460e-01     0.3848\n",
       "        3    -14.56058094      7.7052e-02     0.3856\n",
       "        4    -14.56582672      3.0375e-02     0.3848\n",
       "        5    -14.56662378      1.1794e-02     0.3842\n",
       "        6    -14.56676403      1.7857e-06     0.3837\n",
       "        7    -14.56676404      1.7598e-08     0.3837\n",
       "        8    -14.56676403      1.3644e-09     0.3837\n",
       "    ------------------------------------------------------------\n",
       "\n",
       "\n",
       "    ============================================================\n",
       "               One-electron energy  -19.1171562840\n",
       "                           Kinetic   14.6787428206\n",
       "                         Potential  -33.7958991046\n",
       "      --------------------------------------------------------  \n",
       "               Two-electron energy    4.5503922505\n",
       "                           Coulomb    7.2329463384\n",
       "                          Exchange   -2.6825540879\n",
       "      --------------------------------------------------------  \n",
       "                           HF HOMO   -8.1986650549\n",
       "                           HF LUMO    2.2431802579\n",
       "                      HF LUMO-LUMO   10.4418453128\n",
       "      --------------------------------------------------------  \n",
       "               Electronic   energy  -14.5667640335\n",
       "               Nuclear   repulsion    0.0000000000\n",
       "               Hartree-Fock energy  -14.5667640335\n",
       "    ============================================================\n",
       "    \n",
       "\n"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val mo_basis : MOBasis.t =\n",
       "  Eigenvalues:    -4.706891    -0.301295     0.082435     0.082435     0.082435 \n",
       "            -- 1 --     -- 2 --    -- 3 --   -- 4 --      -- 5 --\n",
       "      1    0.997996    -0.22266          0         0           -0\n",
       "      2   0.0161427    0.288767 2.2404E-15         0 -7.73804E-16\n",
       "      3           0 3.58344E-16   0.258415  0.037609    0.0381979\n",
       "                ...         ...        ...       ...          ...\n",
       "      7           0           0    0.78727  0.114442     0.116234\n",
       "      8           0           0  -0.128237 0.0798944     0.789475\n",
       "      9           0           0  -0.100874  0.791587   -0.0964853\n",
       "  \n",
       "  Eigenvalues:     0.439754     0.464931     0.464931     0.464931 \n",
       "              -- 6 --      -- 7 --     -- 8 --      -- 9 --\n",
       "      1    0.00338551            0          -0            0\n",
       "      2       2.01724 -3.80748E-15 4.70544E-15 -6.99768E-16\n",
       "      3   1.56069E-15        1.292    0.249105    0.0981672\n",
       "                  ...          ...         ...          ...\n",
       "      7   -2.3622E-15     -1.05652    -0.20374   -0.0802898\n",
       "      8   1.64398E-15     0.172096   -0.528064     -0.92522\n",
       "      9   2.61906E-15     0.135375   -0.918775     0.549573\n",
       "  \n",
       "  \n"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "let hf =  HartreeFock.make ~guess:`Hcore  simulation  ;;\n",
    "\n",
    "let mo_basis = MOBasis.of_hartree_fock hf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FCI-F12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notations:\n",
    "\n",
    "* $\\langle ij || kl \\rangle = \\int \\phi_i(r_1) \\phi_j(r_2) \\frac{1}{r_{12}} \\phi_k(r1) \\phi_l(r2) $ \n",
    "* $\\left[ ij || kl \\right] = \\int \\phi_i(r_1) \\phi_j(r_2) f_{12} \\phi_k(r1) \\phi_l(r2) $ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Common functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "val f12 : F12factor.t =\n",
       "  {F12factor.expo_s = 1.;\n",
       "   gaussian =\n",
       "    {GaussianOperator.coef_g =\n",
       "      [  -0.314400  -0.303700  -0.168100  -0.098110  -0.060240  -0.037260 ];\n",
       "     expo_sg =\n",
       "      [   0.220900   1.004000   3.622000  12.160000  45.870000 254.400000 ];\n",
       "     expo_sg_inv =\n",
       "      [   4.526935   0.996016   0.276091   0.082237   0.021801   0.003931 ]}}\n"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val mo_num : int = 9\n"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val pp_spindet : Format.formatter -> Spindeterminant.t -> unit = <fun>\n"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val pp_det : Format.formatter -> Determinant.t -> unit = <fun>\n"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "let f12 = Util.of_some @@ Simulation.f12 simulation \n",
    "\n",
    "let mo_num = MOBasis.size mo_basis \n",
    "\n",
    "let pp_spindet = Spindeterminant.pp mo_num\n",
    "\n",
    "let pp_det = Determinant.pp mo_num\n",
    "\n",
    ";;\n",
    "\n",
    "#install_printer pp_spindet ;;\n",
    "#install_printer pp_det ;;\n",
    "\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "let simulation_aux = \n",
    "    let charge        = Charge.to_int @@ Simulation.charge simulation \n",
    "    and multiplicity  = Electrons.multiplicity @@ Simulation.electrons simulation\n",
    "    and nuclei        = Simulation.nuclei simulation\n",
    "    in\n",
    "    let general_basis = \n",
    "      Basis.general_basis @@ Simulation.basis simulation\n",
    "    in\n",
    "    GeneralBasis.combine [\n",
    "      general_basis ; GeneralBasis.read aux_basis_filename\n",
    "    ]\n",
    "    |> Basis.of_nuclei_and_general_basis nuclei\n",
    "    |> Simulation.make ~f12 ~charge ~multiplicity ~nuclei \n",
    "\n",
    "\n",
    "let aux_basis = \n",
    "  MOBasis.of_mo_basis simulation_aux mo_basis\n",
    "  \n",
    "let () = ignore @@ MOBasis.f12_ints   aux_basis\n",
    "let () = ignore @@ MOBasis.two_e_ints aux_basis\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determinant-based functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\langle I | \\hat{H} | J \\rangle = \\begin{cases}\n",
    "\\sum_i h_{ii} + \\frac{1}{2} \\sum_{ij} \\langle ij || ij \\rangle \\text{ when } |J\\rangle = |I\\rangle \\\\\n",
    "h_{ik} + \\sum_{j} \\langle ij || kj \\rangle \\text{ when } |J\\rangle = \\hat{T}_i^k |I\\rangle \\\\\n",
    "\\langle ij || kl \\rangle \\text{ when } |J\\rangle = \\hat{T}_{ij}^{kl} |I\\rangle \\\\\n",
    "\\end{cases}$\n",
    "\n",
    "\n",
    "$\\langle I | \\hat{F} | J \\rangle = \\begin{cases}\n",
    "\\sum_i f_{ii} + \\frac{1}{2} \\sum_{ij} \\langle ij || ij \\rangle \\text{ when } |J\\rangle = |I\\rangle \\\\\n",
    "f_{ik} + \\sum_{j} \\langle ij || kj \\rangle \\text{ when } |J\\rangle = \\hat{T}_i^k |I\\rangle \\\\\n",
    "\\langle ij || kl \\rangle \\text{ when } |J\\rangle = \\hat{T}_{ij}^{kl} |I\\rangle \\\\\n",
    "\\end{cases}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "val cancel_singles : bool = false\n"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val f12_integrals :\n",
       "  MOBasis.t ->\n",
       "  ('a -> 'b -> 'c -> float) *\n",
       "  (int -> int -> int -> int -> 'd -> 'd -> float) * 'e option = <fun>\n"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val h_ij : MOBasis.t -> CIMatrixElement.De.t -> CIMatrixElement.De.t -> float =\n",
       "  <fun>\n"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val f_ij : MOBasis.t -> CIMatrixElement.De.t -> CIMatrixElement.De.t -> float =\n",
       "  <fun>\n"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val hf_ij :\n",
       "  MOBasis.t -> CIMatrixElement.De.t -> CIMatrixElement.De.t -> float list =\n",
       "  <fun>\n"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "let cancel_singles = false \n",
    "\n",
    "let f12_integrals mo_basis =\n",
    "\n",
    "  \n",
    "  let two_e_ints  = MOBasis.f12_ints mo_basis in\n",
    "  let f2 i j k l s s' =\n",
    "      let ijkl = F12.get_phys two_e_ints i j k l in\n",
    "      let ijlk = F12.get_phys two_e_ints i j l k in\n",
    "      if s' <> s then\n",
    "        ijkl \n",
    "      else\n",
    "        (ijkl -. ijlk) \n",
    "  in\n",
    "  ( (fun _ _ _ -> 0.),\n",
    "    (if cancel_singles then\n",
    "        (fun i j k l s s' ->\n",
    "        (* Required to cancel out single excitations *)\n",
    "           if (i=k && j<>l) || (j=l && i<>k) then\n",
    "              0.\n",
    "           else\n",
    "              f2 i j k l s s'\n",
    "        ) \n",
    "    else f2),\n",
    "    None\n",
    "  )\n",
    "\n",
    "let h_ij mo_basis ki kj =\n",
    "  let integrals =\n",
    "    List.map (fun f -> f mo_basis)\n",
    "      [ CI.h_integrals ]\n",
    "  in\n",
    "  CIMatrixElement.make integrals ki kj \n",
    "  |> List.hd\n",
    "\n",
    "\n",
    "let f_ij mo_basis ki kj =\n",
    "  let integrals =\n",
    "    List.map (fun f -> f mo_basis)\n",
    "      [ f12_integrals ]\n",
    "  in\n",
    "  CIMatrixElement.make integrals ki kj \n",
    "  |> List.hd                                                                                  \n",
    "\n",
    "\n",
    "let hf_ij mo_basis ki kj =\n",
    "  let integrals =\n",
    "    List.map (fun f -> f mo_basis)\n",
    "      [ CI.h_integrals ; f12_integrals ]\n",
    "  in\n",
    "  CIMatrixElement.make integrals ki kj\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determinant space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "- : Spindeterminant.t =\n",
       "+1 +-+-+-+---------------------------------------------------------\n"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Spindeterminant.of_list 10 [1; 3 ; 5; 7]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "val det_space : SpindeterminantSpace.t =\n",
       "   [       0 +1 ++--------------------------------------------------------------\n",
       "           1 +1 --+-------------------------------------------------------------\n",
       "           2 +1 ---+------------------------------------------------------------\n",
       "           3 +1 ----+-----------------------------------------------------------\n",
       "           4 +1 +----+----------------------------------------------------------\n",
       "           5 +1 -+---+----------------------------------------------------------\n",
       "           6 +1 --+--+----------------------------------------------------------\n",
       "           7 +1 ---+-+----------------------------------------------------------\n",
       "           8 +1 ----++----------------------------------------------------------\n",
       "           9 +1 ------+---------------------------------------------------------\n",
       "          10 +1 +++----+--------------------------------------------------------\n",
       "          11 +1 ---+---+--------------------------------------------------------\n",
       "          12 +1 ----+--+--------------------------------------------------------\n",
       "          13 +1 -----+-+--------------------------------------------------------\n",
       "          14 +1 ++----++--------------------------------------------------------\n",
       "          15 +1 --+---++--------------------------------------------------------\n",
       "          16 +1 ---+--++--------------------------------------------------------\n",
       "          17 +1 ----+-++--------------------------------------------------------\n",
       "          18 +1 +----+++--------------------------------------------------------\n",
       "          19 +1 -+---+++--------------------------------------------------------\n",
       "          20 +1 --+--+++--------------------------------------------------------\n",
       "          21 +1 ---+-+++--------------------------------------------------------\n",
       "          22 +1 ----++++--------------------------------------------------------\n",
       "          23 +1 +++-----+-------------------------------------------------------\n",
       "          24 +1 ---+----+-------------------------------------------------------\n",
       "          25 +1 ----+---+-------------------------------------------------------\n",
       "          26 +1 +----+--+-------------------------------------------------------\n",
       "          27 +1 -+---+--+-------------------------------------------------------\n",
       "          28 +1 --+--+--+-------------------------------------------------------\n",
       "          29 +1 ---+-+--+-------------------------------------------------------\n",
       "          30 +1 ----++--+-------------------------------------------------------\n",
       "          31 +1 +-----+-+-------------------------------------------------------\n",
       "          32 +1 -+----+-+-------------------------------------------------------\n",
       "          33 +1 --+---+-+-------------------------------------------------------\n",
       "          34 +1 ---+--+-+-------------------------------------------------------\n",
       "          35 +1 ----+-+-+-------------------------------------------------------\n",
       "    ]\n"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "let det_space = \n",
    "(*\n",
    "    DeterminantSpace.fci_f12_of_mo_basis  aux_basis  ~frozen_core  mo_num\n",
    "    *)\n",
    "    SpindeterminantSpace.fci_of_mo_basis  mo_basis  ~frozen_core (Electrons.n_alfa @@ Simulation.electrons simulation)\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "unsigned int v; // current permutation of bits \n",
    "unsigned int w; // next permutation of bits\n",
    "\n",
    "unsigned int t = v | (v - 1); // t gets v's least significant 0 bits set to 1\n",
    "// Next set to 1 the most significant bit to change, \n",
    "// set to 0 the least significant ones, and add the necessary 1 bits.\n",
    "w = (t + 1) | (((~t & -~t) - 1) >> (__builtin_ctz(v) + 1));  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "- : int = 31\n"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "let _ = \n",
    "    let open Bitstring in\n",
    "(*\n",
    "\n",
    "    let u = Bitstring.of_int 63 in\n",
    "\n",
    "    let t  = logor u @@ minus_one u in\n",
    "    let not_t = lognot t in\n",
    "    let neg_not_t = neg not_t in\n",
    "    \n",
    "    shift_right\n",
    "        (minus_one @@ logand not_t neg_not_t) \n",
    "        ( (trailing_zeros u) + 1),\n",
    "        *)\n",
    "    popcount (Bitstring.of_int 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "val elec_num : int = 4\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val alfa_num : int = 2\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val mo_class : MOClass.t =\n",
       "  [Active 1;Active 2;Active 3;Active 4;Active 5;Active 6;Active 7;Active 8;\n",
       "  Active 9]\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val m : int list -> Bitstring.t = <fun>\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val occ_mask : Bitstring.t =\n",
       "  ----------------------------------------------------------------\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val active_mask : Bitstring.t =\n",
       "  +++++++++-------------------------------------------------------\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val neg_active_mask : Bitstring.t =\n",
       "  ---------+++++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val permtutations : int -> int -> Bitstring.t list = <fun>\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "val spin_determinants : Bitstring.t list =\n",
       "  [++--------------------------------------------------------------;\n",
       "   --+-------------------------------------------------------------;\n",
       "   ---+------------------------------------------------------------;\n",
       "   ----+-----------------------------------------------------------;\n",
       "   -----+----------------------------------------------------------;\n",
       "   ++----+---------------------------------------------------------;\n",
       "   --+---+---------------------------------------------------------;\n",
       "   ---+--+---------------------------------------------------------;\n",
       "   ----+-+---------------------------------------------------------;\n",
       "   +----++---------------------------------------------------------;\n",
       "   -+---++---------------------------------------------------------;\n",
       "   --+--++---------------------------------------------------------;\n",
       "   ---+-++---------------------------------------------------------;\n",
       "   ----+++---------------------------------------------------------;\n",
       "   +------+--------------------------------------------------------;\n",
       "   -+-----+--------------------------------------------------------;\n",
       "   --+----+--------------------------------------------------------;\n",
       "   ---+---+--------------------------------------------------------;\n",
       "   ----+--+--------------------------------------------------------;\n",
       "   +----+-+--------------------------------------------------------;\n",
       "   -+---+-+--------------------------------------------------------;\n",
       "   --+--+-+--------------------------------------------------------;\n",
       "   ---+-+-+--------------------------------------------------------;\n",
       "   ----++-+--------------------------------------------------------;\n",
       "   +-----++--------------------------------------------------------;\n",
       "   -+----++--------------------------------------------------------;\n",
       "   --+---++--------------------------------------------------------;\n",
       "   ---+--++--------------------------------------------------------;\n",
       "   ----+-++--------------------------------------------------------;\n",
       "   +----+++--------------------------------------------------------;\n",
       "   -+---+++--------------------------------------------------------;\n",
       "   --+--+++--------------------------------------------------------;\n",
       "   ---+-+++--------------------------------------------------------;\n",
       "   ----++++--------------------------------------------------------;\n",
       "   +-------+-------------------------------------------------------;\n",
       "   -+------+-------------------------------------------------------]\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "let elec_num = Electrons.n_elec @@ Simulation.electrons simulation\n",
    "\n",
    "let alfa_num = Electrons.n_alfa @@ Simulation.electrons simulation\n",
    "\n",
    "let mo_class = MOClass.fci ~frozen_core mo_basis \n",
    "\n",
    "\n",
    "let m l = \n",
    "    List.fold_left (fun accu i -> let j = i-1 in\n",
    "      Bitstring.logor accu (Bitstring.shift_left_one mo_num j)\n",
    "                   ) (Bitstring.zero mo_num) l\n",
    "                   \n",
    "let occ_mask    = m (MOClass.core_mos    mo_class)\n",
    "\n",
    "let active_mask = m (MOClass.active_mos  mo_class)\n",
    "\n",
    "\n",
    "let neg_active_mask = Bitstring.lognot active_mask \n",
    "\n",
    "  (* Here we generate the FCI space and filter out unwanted determinants\n",
    "     with excitations involving the core electrons. This should be improved. *)\n",
    "     \n",
    "let permtutations m n =                                                                   \n",
    "  let open Bitstring in\n",
    "  let rec aux k u rest =\n",
    "    if k=1 then\n",
    "      List.rev (u :: rest)\n",
    "    else\n",
    "      let t   = (logor u (minus_one u)) in\n",
    "      let t'  = plus_one t in\n",
    "      let t'' = shift_right (minus_one (logand (lognot t) (neg @@ lognot t))) (trailing_zeros u + 1) in                                                             \n",
    "      (*\n",
    "      let t'' = shift_right (minus_one (logand (lognot t) t')) (trailing_zeros u + 1) in\n",
    "      *)\n",
    "      (aux [@tailcall]) (k-1) (logor t' t'') (u :: rest)\n",
    "  in\n",
    "  aux (Util.binom n m) (minus_one (shift_left_one n m)) []\n",
    "\n",
    "     \n",
    "let spin_determinants =\n",
    "    Bitstring.permtutations alfa_num mo_num\n",
    "    (*\n",
    "    |> List.filter (fun b -> Bitstring.logand neg_active_mask b = occ_mask)\n",
    "    |> List.map (fun b -> Spindeterminant.of_bitstring b)\n",
    "    |> Array.of_list\n",
    "    *)\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "let ci = CI.make ~n_states:state det_space\n",
    "\n",
    "let ci_coef, ci_energy = Lazy.force ci.eigensystem "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creates a function `f : Determinant.t -> bool` that returns `true` when the determinant has one or two electrons in the CABS."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "let is_a_double det_space =\n",
    "    let mo_class      = DeterminantSpace.mo_class det_space in\n",
    "    let mo_num = Array.length @@ MOClass.mo_class_array mo_class in\n",
    "    let m l =\n",
    "      List.fold_left (fun accu i ->\n",
    "          let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)\n",
    "        ) (Bitstring.zero mo_num) l\n",
    "    in\n",
    "    let aux_mask   = m (MOClass.auxiliary_mos mo_class) in\n",
    "    fun k ->                                                                                  \n",
    "      let alfa =\n",
    "        Determinant.alfa k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let beta =\n",
    "        Determinant.beta k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let a = Bitstring.logand aux_mask alfa\n",
    "      and b = Bitstring.logand aux_mask beta\n",
    "      in\n",
    "      match Bitstring.popcount a + Bitstring.popcount b with\n",
    "      | 2 | 1 -> true\n",
    "      | 0 | _ -> false\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Permutation operator $p_{12}$ that generates a new determinant with electrons 1 and 2 swapped."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "let p12 det_space =\n",
    "    let mo_class      = DeterminantSpace.mo_class det_space in\n",
    "    let mo_num = Array.length @@ MOClass.mo_class_array mo_class in\n",
    "    let m l =\n",
    "      List.fold_left (fun accu i ->\n",
    "          let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)\n",
    "        ) (Bitstring.zero mo_num) l\n",
    "    in\n",
    "    let aux_mask = m (MOClass.auxiliary_mos mo_class) in\n",
    "    let not_aux_mask =\n",
    "      Bitstring.(shift_left_one mo_num (mo_num-1) |> minus_one |> logxor aux_mask)\n",
    "    in\n",
    "    fun k ->\n",
    "      let alfa =\n",
    "        Determinant.alfa k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let beta =\n",
    "        Determinant.beta k\n",
    "        |> Spindeterminant.bitstring\n",
    "      in\n",
    "      let a = Bitstring.logand aux_mask alfa\n",
    "      and b = Bitstring.logand aux_mask beta\n",
    "      in\n",
    "      match Bitstring.popcount a, Bitstring.popcount b with\n",
    "      | 2, 0  \n",
    "      | 0, 2 -> Some (Determinant.negate_phase k)\n",
    "      | 1, 1 -> Some (Determinant.of_spindeterminants\n",
    "                  (Spindeterminant.of_bitstring  @@\n",
    "                    Bitstring.(logor b (logand not_aux_mask alfa)) )\n",
    "                  (Spindeterminant.of_bitstring  @@\n",
    "                    Bitstring.(logor a (logand not_aux_mask beta))\n",
    "                  ) )\n",
    "      | 1, 0 \n",
    "      | 0, 1 -> Some k\n",
    "      | _ -> None\n",
    "                    \n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determinants of the FCI space"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "let in_dets =     \n",
    "    DeterminantSpace.determinant_stream ci.CI.det_space \n",
    "    |> Util.stream_to_list\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integral-based functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "OCaml 4.07.1",
   "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.07.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}