diff --git a/.ocamlinit b/.ocamlinit
index 7d69cc6..af6a68b 100644
--- a/.ocamlinit
+++ b/.ocamlinit
@@ -5,12 +5,21 @@ let () =
 
 #use "topfind";;
 #require "lacaml";;
+#require "alcotest";;
+#require "zarith";;
+#require "getopt";;
 #directory "_build";;
-#directory "_build/Utils";;
-#directory "_build/Nuclei";;
 #directory "_build/Basis";;
-#directory "_build/HartreeFock";;
+#directory "_build/CI";;
+#directory "_build/MOBasis";;
+#directory "_build/Nuclei";;
+#directory "_build/Parallel";;
+#directory "_build/Perturbation";;
+#directory "_build/SCF";;
+#directory "_build/Utils";;
+
+#load "Constants.cmo";;
+#load_rec "Util.cma";;
 
-#use "topfind";;
 
 
diff --git a/Makefile.include b/Makefile.include
index 7eadff4..3747291 100644
--- a/Makefile.include
+++ b/Makefile.include
@@ -1,9 +1,8 @@
 .NOPARALLEL:
 
-INCLUDE_DIRS=Parallel,Nuclei,Utils,Basis,SCF,MOBasis,CI,F12,Perturbation
 LIBS=
 PKGS=
-OCAMLBUILD=ocamlbuild -j 0 -cflags $(ocamlcflags) -lflags $(ocamllflags) $(ocamldocflags) -Is $(INCLUDE_DIRS) -ocamlopt $(ocamloptflags) $(mpi) 
+OCAMLBUILD=ocamlbuild -j 0 -cflags $(ocamlcflags) -lflags $(ocamllflags) $(ocamldocflags) -ocamlopt $(ocamloptflags) $(mpi) 
 
 MLLFILES=$(filter-out $(wildcard _build/*), $(wildcard */*.mll) $(wildcard *.mll)) Utils/math_functions.c
 MLYFILES=$(filter-out $(wildcard _build/*), $(wildcard */*.mly) $(wildcard *.mly))
@@ -12,7 +11,7 @@ MLIFILES=$(filter-out $(wildcard Parallel_*/*) $(wildcard _build/*), $(wildcard
 
 ALL_NATIVE=$(patsubst %.ml,%.native,$(wildcard run_*.ml)) 
 ALL_BYTE=$(patsubst %.ml,%.byte,$(wildcard run_*.ml)) 
-ALL_EXE=$(ALL_BYTE) $(ALL_NATIVE)
+ALL_EXE=$(ALL_BYTE) $(ALL_NATIVE) lib
 
 .PHONY: default doc
 
@@ -21,8 +20,10 @@ default: $(ALL_EXE)
 
 tests: run_tests.native
 
-bytelib:
-	ocamlmklib -o bytelib ./_build/Utils/Util.cmo 
+lib: _build/Utils/Util.cma
+
+_build/Utils/Util.cma: _build/Utils/Util.cmo _build/Utils/math_functions.o
+	ocamlmklib -o _build/Utils/Util _build/Utils/Util.cmo _build/Utils/math_functions.o 
 
 QCaml.odocl: $(MLIFILES)
 	ls $(MLIFILES) | sed "s/\.mli//" > QCaml.odocl
diff --git a/Notebooks/F12_matrix.ipynb b/Notebooks/F12_matrix.ipynb
new file mode 100644
index 0000000..4dee672
--- /dev/null
+++ b/Notebooks/F12_matrix.ipynb
@@ -0,0 +1,159 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Test of F12 matrix elements"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "- : unit = ()\n",
+      "Findlib has been successfully loaded. Additional directives:\n",
+      "  #require \"package\";;      to load a package\n",
+      "  #list;;                   to list the available packages\n",
+      "  #camlp4o;;                to load camlp4 (standard syntax)\n",
+      "  #camlp4r;;                to load camlp4 (revised syntax)\n",
+      "  #predicates \"p,q,...\";;   to set these predicates\n",
+      "  Topfind.reset();;         to force that packages will be reloaded\n",
+      "  #thread;;                 to enable threads\n",
+      "\n",
+      "- : unit = ()\n"
+     ]
+    }
+   ],
+   "source": [
+    "#cd \"/home/scemama/QCaml\";;\n",
+    "#use \"topfind\";;\n",
+    "#require \"jupyter.notebook\";;\n",
+    "\n",
+    "#require \"lacaml\";;\n",
+    "#require \"alcotest\";;\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",
+    "#load \"Constants.cmo\";;\n",
+    "#load_rec \"Util.cma\";;\n",
+    "#load_rec \"Matrix.cmo\";;\n",
+    "#load_rec \"Simulation.cmo\";;\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Simulation\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "File \"[12]\", line 3, characters 4-22:\n",
+      "Warning 26: unused variable aux_basis_filename.\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "- : Simulation.t = <abstr>\n"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "let basis_file  = \"/home/scemama/qp2/data/basis/6-31g\" in\n",
+    "          \n",
+    "let aux_basis_filename = \"/home/scemama/qp2/data/basis/cc-pvdz\" in\n",
+    "  \n",
+    "let nuclei_file = \"be.xyz\" in\n",
+    "  \n",
+    "let charge = 0 in\n",
+    "\n",
+    "let multiplicity = 1 in\n",
+    "\n",
+    "let f12 = F12factor.gaussian_geminal 1.0 in\n",
+    "\n",
+    "let simulation =\n",
+    "  Simulation.of_filenames ~f12 ~charge ~multiplicity ~nuclei:nuclei_file basis_file \n",
+    "in\n",
+    "simulation\n",
+    "  \n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "NameError",
+     "evalue": "name 'version' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
+      "\u001b[0;32m<ipython-input-10-605b5d1778ad>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mversion\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
+      "\u001b[0;31mNameError\u001b[0m: name 'version' is not defined"
+     ]
+    }
+   ],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "OCaml default",
+   "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
+}
diff --git a/Notebooks/SpVec.ipynb b/Notebooks/SpVec.ipynb
index 29df2ae..4070396 100644
--- a/Notebooks/SpVec.ipynb
+++ b/Notebooks/SpVec.ipynb
@@ -74,19 +74,6 @@
     "----------------"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "A sparse vector is a structure made of:\n",
-    "* The dimension of the vector space\n",
-    "* The number of non-zeros\n",
-    "* An array of indices\n",
-    "* An array of values\n",
-    "\n",
-    "The indices are stored in an ``int Bigarray`` and the values are stored in a ``Lacaml.Vec.t``\n"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -147,7 +134,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 32,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
     {
@@ -156,7 +143,7 @@
        "val pp : Format.formatter -> t -> unit = <fun>\n"
       ]
      },
-     "execution_count": 32,
+     "execution_count": 3,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -189,7 +176,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
@@ -198,7 +185,7 @@
        "val make0 : int -> t = <fun>\n"
       ]
      },
-     "execution_count": 5,
+     "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -214,7 +201,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
     {
@@ -223,7 +210,7 @@
        "val of_vec : ?threshold:float -> Lacaml.D.vec -> t = <fun>\n"
       ]
      },
-     "execution_count": 6,
+     "execution_count": 5,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -262,7 +249,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
@@ -271,7 +258,7 @@
        "val of_array : ?threshold:float -> float array -> t = <fun>\n"
       ]
      },
-     "execution_count": 7,
+     "execution_count": 6,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -284,7 +271,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
@@ -293,7 +280,7 @@
        "val copy : t -> t = <fun>\n"
       ]
      },
-     "execution_count": 8,
+     "execution_count": 7,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -312,7 +299,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [
     {
@@ -321,7 +308,7 @@
        "val map : ?threshold:float -> (float -> float) -> t -> t = <fun>\n"
       ]
      },
-     "execution_count": 37,
+     "execution_count": 8,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -361,16 +348,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 38,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "val x : t = {dim: 10 }\n"
+       "val x : t =\n",
+       "  {dim = 10; nnz = 0; indices = <abstr>;\n",
+       "   values =\n",
+       "              R1           R2 R3              R30          R31          R32\n",
+       "    6.91705E-310 6.91705E-310  0 ... 6.91705E-310 6.91705E-310 6.91705E-310}\n"
       ]
      },
-     "execution_count": 38,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -381,17 +372,20 @@
        "                              1 -2  0 ...  0  3  0\n"
       ]
      },
-     "execution_count": 38,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     },
     {
      "data": {
       "text/plain": [
-       "val sparse_a : t = {dim: 9 (1, 1) (2, -2) (5, 0.5) (6, 1e-08) (8, 3) }\n"
+       "val sparse_a : t =\n",
+       "  {dim = 9; nnz = 5; indices = <abstr>;\n",
+       "   values = R1 R2  R3    R4 R5\n",
+       "             1 -2 0.5 1E-08  3}\n"
       ]
      },
-     "execution_count": 38,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -401,7 +395,7 @@
        "- : bool = true\n"
       ]
      },
-     "execution_count": 38,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -411,7 +405,7 @@
        "- : bool = false\n"
       ]
      },
-     "execution_count": 38,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -440,7 +434,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 39,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
@@ -449,7 +443,7 @@
        "val dim : t -> int = <fun>\n"
       ]
      },
-     "execution_count": 39,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -459,7 +453,7 @@
        "val nnz : t -> int = <fun>\n"
       ]
      },
-     "execution_count": 39,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -471,7 +465,7 @@
        "  <fun>\n"
       ]
      },
-     "execution_count": 39,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -481,7 +475,7 @@
        "val values : t -> L.Vec.t = <fun>\n"
       ]
      },
-     "execution_count": 39,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -491,7 +485,7 @@
        "val density : t -> float = <fun>\n"
       ]
      },
-     "execution_count": 39,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -506,7 +500,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 40,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [
     {
@@ -515,7 +509,7 @@
        "val get : t -> int -> float = <fun>\n"
       ]
      },
-     "execution_count": 40,
+     "execution_count": 11,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -549,7 +543,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 41,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [
     {
@@ -558,7 +552,7 @@
        "val iter : (int -> float -> 'a) -> t -> unit = <fun>\n"
       ]
      },
-     "execution_count": 41,
+     "execution_count": 12,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -579,7 +573,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 42,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [
     {
@@ -588,7 +582,7 @@
        "- : bool = true\n"
       ]
      },
-     "execution_count": 42,
+     "execution_count": 13,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -609,7 +603,7 @@
        "- : unit = ()\n"
       ]
      },
-     "execution_count": 42,
+     "execution_count": 13,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -628,7 +622,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 43,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [
     {
@@ -637,7 +631,7 @@
        "val to_assoc_list : t -> (int * float) list = <fun>\n"
       ]
      },
-     "execution_count": 43,
+     "execution_count": 14,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -647,7 +641,7 @@
        "val to_vec : t -> Lacaml.D.vec = <fun>\n"
       ]
      },
-     "execution_count": 43,
+     "execution_count": 14,
      "metadata": {},
      "output_type": "execute_result"
     }
diff --git a/_tags b/_tags
index 0f8a66a..94195d5 100644
--- a/_tags
+++ b/_tags
@@ -1,4 +1,13 @@
-true: package(str,unix,bigarray,lacaml,alcotest,zarith,getopt)
-<*.byte>  : linkdep(Utils/math_functions.o), custom
-<*.native>: linkdep(Utils/math_functions.o)
-<odoc-ltxhtml>: not_hygienic
+not <**/*.cma>   : linkdep(Utils/math_functions.o),  package(str,unix,bigarray,lacaml,alcotest,zarith,getopt)
+<**/*.byte>      : custom
+<**/*.native>    : linkdep(Utils/math_functions.o)
+<odoc-ltxhtml>   : not_hygienic
+<Parallel>       : include
+<Nuclei>         : include
+<Utils>          : include
+<Basis>          : include
+<SCF>            : include
+<MOBasis>        : include
+<CI>             : include
+<F12>            : include
+<Perturbation>   : include
diff --git a/configure b/configure
index 905185b..bdbac95 100755
--- a/configure
+++ b/configure
@@ -11,9 +11,9 @@ doc='$(prefix)/doc'
 share='$(prefix)/share'
 man='$(prefix)/man'
 etc='$(prefix)/etc'
-ocamlcflags='"-g -warn-error A -cc gcc -ccopt -Ofast -ccopt -march=native "'
+ocamlcflags='"-g -warn-error A -cc gcc -ccopt -Ofast -ccopt -march=native -ccopt=-fPIC "'
 ocamllflags='"-g -warn-error A"'
-ocamloptflags='"opt -O3 -nodynlink -remove-unused-arguments -rounds 16 -inline 100 -inline-max-unroll 100 -cc gcc -ccopt -Ofast -ccopt -march=native"'
+ocamloptflags='"opt -O3 -nodynlink -remove-unused-arguments -rounds 16 -inline 100 -inline-max-unroll 100 -cc gcc -ccopt -Ofast -ccopt -march=native -ccopt=-fPIC "'
 ocamldocflags='-docflags "-g ltxhtml.cma -sort -css-style $(PWD)/style.css -colorize-code"'
 mpi=''
 ocamloptflags_profile='"opt -p -cc gcc -ccopt -Ofast -ccopt -march=native -ccopt -no-pie "'