diff --git a/ao_basis/lib/ao_basis.ml b/ao_basis/lib/ao_basis.ml
new file mode 100644
index 0000000..f3cf2e3
--- /dev/null
+++ b/ao_basis/lib/ao_basis.ml
@@ -0,0 +1,8 @@
+type basis =
+  | Unknown
+  | Gaussian of Ao_basis_gaussian.t
+
+type t =
+  { basis : basis  ;
+    cartesian : bool
+  }
diff --git a/ao_basis/lib/ao_basis.mli b/ao_basis/lib/ao_basis.mli
new file mode 100644
index 0000000..a1a7dd9
--- /dev/null
+++ b/ao_basis/lib/ao_basis.mli
@@ -0,0 +1,10 @@
+(** Data structure for Atomic Orbitals. *)
+
+type basis =
+  | Unknown
+  | Gaussian of Ao_basis_gaussian.t
+
+type t =
+  { basis : basis  ;
+    cartesian : bool
+  }
diff --git a/ao_basis/lib/ao_basis_gaussian.ml b/ao_basis/lib/ao_basis_gaussian.ml
new file mode 100644
index 0000000..4834cf8
--- /dev/null
+++ b/ao_basis/lib/ao_basis_gaussian.ml
@@ -0,0 +1,189 @@
+module GB = Qcaml_gaussian_basis
+module GI = Qcaml_gaussian_integrals
+  
+type t =
+{
+  basis              :  GB.Basis.t         ;
+  overlap            :  GI.Overlap.t             lazy_t;
+  multipole          :  GI.Multipole.t           lazy_t;
+  ortho              :  GI.Orthonormalization.t  lazy_t;
+  eN_ints            :  GI.Electron_nucleus.t    lazy_t;
+  kin_ints           :  GI.Kinetic.t             lazy_t;
+  ee_ints            :  GI.Eri.t                 lazy_t;
+  ee_lr_ints         :  GI.Eri_long_range.t      lazy_t;
+  f12_ints           :  GI.F12.t                 lazy_t;
+  f12_over_r12_ints  :  GI.Screened_eri.t        lazy_t;
+}
+
+(*
+let basis        t = t.basis
+let overlap      t = Lazy.force t.overlap
+let multipole    t = Lazy.force t.multipole
+let ortho        t = Lazy.force t.ortho
+let eN_ints      t = Lazy.force t.eN_ints
+let kin_ints     t = Lazy.force t.kin_ints
+let ee_ints      t = Lazy.force t.ee_ints
+let ee_lr_ints   t = Lazy.force t.ee_lr_ints
+let f12_ints     t = Lazy.force t.f12_ints
+let f12_over_r12_ints t = Lazy.force t.f12_over_r12_ints
+let cartesian    t = t.cartesian
+
+module Cs = ContractedShell
+
+let values t point =
+  let result = Vec.create (Basis.size t.basis) in
+  Array.iter (fun shell ->
+      Cs.values shell point
+      |> Array.iteri
+           (fun i_c value ->
+             let i = Cs.index shell + i_c + 1 in
+             result.{i} <- value)
+    ) (Basis.contracted_shells t.basis);
+  result
+
+  
+let make ~cartesian ~basis ?f12 nuclei =
+  
+  let overlap =
+    lazy (
+        Overlap.of_basis basis
+      ) in
+  
+  let ortho =
+    lazy (
+        Orthonormalization.make ~cartesian ~basis (Lazy.force overlap)
+      ) in
+
+  let eN_ints  =
+    lazy (
+      Electron_nucleus.of_basis_nuclei ~basis nuclei
+      ) in
+
+  let kin_ints =
+    lazy (
+        Kinetic.of_basis basis
+      ) in
+
+  let ee_ints  =
+    lazy (
+        ERI.of_basis basis
+      ) in
+
+  let ee_lr_ints  =
+    lazy (
+        ERI_lr.of_basis basis
+      ) in
+
+  let f12_ints =
+    lazy (
+        F12.of_basis basis
+      ) in
+
+  let f12_over_r12_ints =
+    lazy (
+        ScreenedERI.of_basis basis
+      ) in
+
+  let multipole =
+    lazy (
+        Multipole.of_basis basis
+      ) in
+
+  { basis ; overlap ; multipole ; ortho ; eN_ints ; kin_ints ; ee_ints ;
+    ee_lr_ints ; f12_ints ; f12_over_r12_ints ; cartesian ;
+  }
+
+
+
+let test_case name t =
+
+  let check_matrix title a r =
+    let a = Mat.to_array a in
+    Mat.to_array r
+    |> Array.iteri (fun i x ->
+          let message =
+            Printf.sprintf "%s line %d" title i
+          in
+          Alcotest.(check (array (float 1.e-10))) message a.(i) x
+        )
+  in
+
+  let check_eri a r =
+    let f { ERI.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } =
+      (i_r1, (j_r2, (k_r1, (l_r2, value))))
+    in
+    let a = ERI.to_list a |> List.rev_map f |> List.rev in
+    let r = ERI.to_list r |> List.rev_map f |> List.rev in
+    Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-10))))))) "ERI" a r
+  in
+
+  let check_eri_lr a r =
+    let f { ERI_lr.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } =
+      (i_r1, (j_r2, (k_r1, (l_r2, value))))
+    in
+    let a = ERI_lr.to_list a |> List.rev_map f |> List.rev in
+    let r = ERI_lr.to_list r |> List.rev_map f |> List.rev in
+    Alcotest.(check (list (pair int (pair int (pair int (pair int (float 1.e-10))))))) "ERI_lr" a r
+  in
+
+  let test_overlap () =
+    let reference =
+      sym_matrix_of_file ("test_files/"^name^"_overlap.ref")
+    in
+    let overlap =
+      Lazy.force t.overlap |> Overlap.matrix
+    in
+    check_matrix "Overlap" overlap reference
+  in
+
+  let test_eN_ints () =
+    let reference =
+      sym_matrix_of_file ("test_files/"^name^"_nuc.ref")
+    in
+    let eN_ints =
+      Lazy.force t.eN_ints |> NucInt.matrix
+    in
+    check_matrix "eN_ints" eN_ints reference
+  in
+
+  let test_kin_ints () =
+    let reference =
+      sym_matrix_of_file ("test_files/"^name^"_kin.ref")
+    in
+    let kin_ints =
+      Lazy.force t.kin_ints |> KinInt.matrix
+    in
+    check_matrix "kin_ints" kin_ints reference
+  in
+
+  let test_ee_ints () =
+    let reference =
+      ERI.of_file ("test_files/"^name^"_eri.ref") ~sparsity:`Dense ~size:(Basis.size t.basis)
+    in
+    let ee_ints =
+      Lazy.force t.ee_ints
+    in
+    check_eri ee_ints reference
+    ;
+  in
+
+  let test_ee_lr_ints () =
+    let reference =
+      ERI_lr.of_file ("test_files/"^name^"_eri_lr.ref") ~sparsity:`Dense
+          ~size:(Basis.size t.basis)
+    in
+    let ee_lr_ints =
+      Lazy.force t.ee_lr_ints
+    in
+    check_eri_lr ee_lr_ints reference
+  in
+
+  [
+    "Overlap",     `Quick,  test_overlap;
+    "eN_ints",     `Quick,  test_eN_ints;
+    "kin_ints",    `Quick,  test_kin_ints;
+    "ee_ints",     `Quick,  test_ee_ints;
+    "ee_lr_ints",  `Quick,  test_ee_lr_ints;
+  ]
+   
+*)
diff --git a/ao_basis/lib/ao_basis_gaussian.mli b/ao_basis/lib/ao_basis_gaussian.mli
new file mode 100644
index 0000000..81e050b
--- /dev/null
+++ b/ao_basis/lib/ao_basis_gaussian.mli
@@ -0,0 +1,69 @@
+(** Data structure for Atomic Orbitals. *)
+
+module GB = Qcaml_gaussian_basis
+module GI = Qcaml_gaussian_integrals
+  
+type t =
+{
+  basis              :  GB.Basis.t         ;
+  overlap            :  GI.Overlap.t             lazy_t;
+  multipole          :  GI.Multipole.t           lazy_t;
+  ortho              :  GI.Orthonormalization.t  lazy_t;
+  eN_ints            :  GI.Electron_nucleus.t    lazy_t;
+  kin_ints           :  GI.Kinetic.t             lazy_t;
+  ee_ints            :  GI.Eri.t                 lazy_t;
+  ee_lr_ints         :  GI.Eri_long_range.t      lazy_t;
+  f12_ints           :  GI.F12.t                 lazy_t;
+  f12_over_r12_ints  :  GI.Screened_eri.t        lazy_t;
+}
+
+
+(*
+(** {1 Accessors} *)
+
+val basis : t -> Basis.t
+(** One-electron basis set *)
+
+val overlap : t -> Overlap.t
+(** Overlap matrix *)
+
+val multipole : t -> Multipole.t
+(** Multipole matrices *)
+
+val ortho : t -> Orthonormalization.t
+(** Orthonormalization matrix of the overlap *)
+
+val eN_ints : t -> NucInt.t
+(** Electron-nucleus potential integrals *)
+
+val ee_ints : t -> ERI.t
+(** Electron-electron potential integrals *)
+
+val ee_lr_ints : t -> ERI_lr.t
+(** Electron-electron long-range potential integrals *)
+
+val f12_ints : t -> F12.t
+(** Electron-electron potential integrals *)
+
+val kin_ints : t -> KinInt.t
+(** Kinetic energy integrals *)
+
+val cartesian : t -> bool
+(** If true, use cartesian Gaussians (6d, 10f, ...) *)
+
+val values : t -> Coordinate.t -> Vec.t
+(** Values of the AOs evaluated at a given point *)
+
+
+
+(** {1 Creators} *)
+
+val make : cartesian:bool -> basis:Basis.t ->  ?f12:F12factor.t ->  Nuclei.t -> t
+(** Creates the data structure for atomic orbitals from a {Basis.t} and the
+    molecular geometry {Nuclei.t} *)
+
+
+(** {2 Tests} *)
+
+val test_case : string -> t -> unit Alcotest.test_case list
+   *)
diff --git a/ao_basis/lib/dune b/ao_basis/lib/dune
new file mode 100644
index 0000000..ae2a316
--- /dev/null
+++ b/ao_basis/lib/dune
@@ -0,0 +1,12 @@
+; name = name of the supermodule that will wrap all source files as submodules
+; public_name = name of the library for ocamlfind and opam
+(library
+ (name qcaml_ao_basis)
+ (public_name qcaml.ao_basis)
+ (libraries
+        qcaml.common
+        qcaml.gaussian_basis
+        qcaml.gaussian_integrals
+	qcaml.operators
+ )
+ (synopsis "Atomic basis set."))
diff --git a/gaussian_integrals/lib/two_electron_integrals.ml b/gaussian_integrals/lib/two_electron_integrals.ml
index 2859a6a..100ba7b 100644
--- a/gaussian_integrals/lib/two_electron_integrals.ml
+++ b/gaussian_integrals/lib/two_electron_integrals.ml
@@ -26,6 +26,7 @@ module type Two_ei_structure =
 module Make(T : Two_ei_structure) = struct
 
   include Four_idx_storage
+  type t = Basis.t Four_idx_storage.t
 
   let class_of_contracted_shell_pair_couple = T.class_of_contracted_shell_pair_couple 
 
diff --git a/gaussian_integrals/lib/two_electron_integrals.mli b/gaussian_integrals/lib/two_electron_integrals.mli
index af36f9e..ba1ed14 100644
--- a/gaussian_integrals/lib/two_electron_integrals.mli
+++ b/gaussian_integrals/lib/two_electron_integrals.mli
@@ -27,8 +27,9 @@ val class_of_contracted_shell_pair_couple :
 module Make : functor (T : Two_ei_structure) -> 
     sig
       include module type of Four_idx_storage
+      type t = Basis.t Four_idx_storage.t
 
-      val of_basis : ?operator:Operator.t ->  Basis.t -> Basis.t t
+      val of_basis : ?operator:Operator.t ->  Basis.t -> t
       (** Compute all ERI's for a given {!Basis.t}. *)
 
     end