diff --git a/Basis/Basis.ml b/Basis/Basis.ml
index d19fe2d..f3488cb 100644
--- a/Basis/Basis.ml
+++ b/Basis/Basis.ml
@@ -1,4 +1,11 @@
-type t = Contracted_shell.t array
+type t = 
+  {
+    size : int;
+    contracted_shells : Contracted_shell.t array;
+  }
+
+let size b = b.size
+let contracted_shells b = b.contracted_shells
 
 (** Returns an array of the basis set per atom *)
 let of_nuclei_and_general_basis n b =
@@ -22,10 +29,19 @@ let of_nuclei_and_general_basis n b =
             (Contracted_shell.index result.(i-1)) +
             (Array.length (Contracted_shell.powers result.(i-1)))) 
     ) result ;
-  result
+
+  let size =
+    let n = ref 0 in
+    for i=0 to (Array.length result) - 1 do
+      n := !n + (Array.length (Contracted_shell.powers (result.(i))))
+    done; !n
+  in
+  { contracted_shells = result ; size }
+
 
 
 let to_string b =
+  let b = b.contracted_shells in
   let line ="
 -----------------------------------------------------------------------
 " in 
@@ -46,21 +62,10 @@ let to_string b =
   )
   ^ line
 
-let file : string option ref = ref None
-
-let set_file f =
-  file := Some f
-
-let general_basis = lazy(
-  match !file with
-  | None -> failwith "basis set file not defined"
-  | Some filename -> Gamess_reader.read ~filename
-)   
-
-
-let basis = lazy (
-  of_nuclei_and_general_basis
-    (Lazy.force Nuclei.nuclei) (Lazy.force general_basis)
-)
+let of_nuclei_and_basis_filename ~nuclei ~filename = 
+  let general_basis = 
+    Gamess_reader.read ~filename
+  in
+  of_nuclei_and_general_basis nuclei general_basis
 
 
diff --git a/Basis/Basis.mli b/Basis/Basis.mli
index 30be9c8..198ed57 100644
--- a/Basis/Basis.mli
+++ b/Basis/Basis.mli
@@ -1,4 +1,10 @@
-type t = Contracted_shell.t array
+type t  
+
+(** Number of contracted Gaussians *)
+val size : t -> int
+
+(** Array of contracted shells *)
+val contracted_shells : t -> Contracted_shell.t array
 
 
 (** Returns an array of the basis set per atom *)
@@ -9,22 +15,7 @@ val of_nuclei_and_general_basis : Nuclei.t -> General_basis.t list -> t
 val to_string : t -> string
 
 
-(** Mutates the state of the file variable to Some f. Required for command-line
-    interface.
-*)
-val set_file : string -> unit
+(** Create a basis using the coordinates of Nuclei and a the filename of
+    the general basis set *)
+val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> filename:string -> t
 
-
-(** Basis set file read and parsed. Requires the set_file function to have
-    been called before.
-*)
-val general_basis :
-  (Element.t * General_basis.general_contracted_shell array) list lazy_t
-
-
-(** Global variable which sets the basis set of the current run.
-    Lazy evaluated from the nuclear coordinates (Nuclei.nuclei) and
-    the basis set file (general_basis). The set_file function has
-    to be called before the basis is used.
-*)
-val basis : Contracted_shell.t array lazy_t
diff --git a/Basis/ERI.ml b/Basis/ERI.ml
index a3623cd..b7324a2 100644
--- a/Basis/ERI.ml
+++ b/Basis/ERI.ml
@@ -4,6 +4,8 @@ open Util
 open Constants
 open Bigarray
 
+type t = (float, float64_elt, fortran_layout) Bigarray.Genarray.t
+
 (** (00|00)^m : Fundamental electron repulsion integral
     $ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $
 
@@ -59,8 +61,7 @@ let index i j k l =
 *)
 
 
-(** Write all integrals to a file with the <ij|kl> convention *)
-let to_file ~filename basis =
+let of_basis basis = 
   let to_int_tuple x = 
     let open Zkey in
     match to_int_tuple Kind_3 x with
@@ -68,19 +69,16 @@ let to_file ~filename basis =
     | _ -> assert false
   in
 
-  let oc = open_out filename in
+  let n =
+    Basis.size basis
+  and shell =
+    Basis.contracted_shells basis
+  in
 
-  Printf.printf "%d shells\n" (Array.length basis);
 
   (* Pre-compute all shell pairs *)
   let shell_pairs =
-    ContractedShellPair.shell_pairs basis
-  in
-  let indices_of_shell_pairs =
-    ContractedShellPair.indices shell_pairs
-  in
-  let unique_shell_pairs =
-    ContractedShellPair.unique shell_pairs
+    ContractedShellPair.shell_pairs shell
   in
 
   (* Pre-compute diagonal integrals for Schwartz *)
@@ -96,8 +94,8 @@ let to_file ~filename basis =
   in
 
   let icount = ref 0 in
-  for i=0 to (Array.length basis) - 1 do
-    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
+  for i=0 to (Array.length shell) - 1 do
+    print_int (Contracted_shell.index shell.(i)) ; print_newline ();
     for j=0 to i do
       let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
         if (schwartz_p_max >= cutoff) then
@@ -108,35 +106,26 @@ let to_file ~filename basis =
 
   (* 4D data initialization *)
   let eri_array = 
-
-    let n = ref 0 in
-    for i=0 to (Array.length basis) - 1 do
-      n := !n + (Array.length (Contracted_shell.powers (basis.(i))))
-    done;
-    let n = !n in
-    Genarray.create Float64 c_layout [| n ; n ; n ; n|]
+    Genarray.create Float64 fortran_layout [| n ; n ; n ; n|]
   in
   Genarray.fill eri_array 0.;
 
   (* Compute ERIs *)
-  let t0 = Unix.gettimeofday () in
 
+  let t0 = Unix.gettimeofday () in
   let inn = ref 0 and out = ref 0 in
 
-  (*
-  for i=0 to (Array.length basis) - 1 do
-    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
+  for i=0 to (Array.length shell) - 1 do
+    print_int (Contracted_shell.index shell.(i)) ; print_newline ();
     for j=0 to i do
-    *)
-  List.iter (fun shell_p ->  
-    let i,j = 
-      Hashtbl.find indices_of_shell_pairs (ContractedShellPair.hash shell_p)
-    in
-    assert (ContractedShellPair.equivalent (ContractedShellPair.create basis.(i) basis.(j)) shell_p);
       let schwartz_p, schwartz_p_max = schwartz.(i).(j) in
       try
         if (schwartz_p_max < cutoff) then raise NullIntegral;
 
+        let
+          shell_p = shell_pairs.(i).(j)
+        in
+
         for k=0 to i do
           for l=0 to k do
             let schwartz_q, schwartz_q_max = schwartz.(k).(l) in
@@ -165,21 +154,18 @@ let to_file ~filename basis =
                   contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q 
             in
 
-      Hashtbl.find_all indices_of_shell_pairs (ContractedShellPair.hash shell_p)
-      |> List.iter (fun (i,j) ->
-            Printf.printf "%d %d\n" i j;
             (* Write the data in the output file *)
             Array.iteri (fun i_c powers_i ->
-                let i_c = (Contracted_shell.index basis.(i)) + i_c  in
+                let i_c = (Contracted_shell.index shell.(i)) + i_c + 1 in
                 let xi = to_int_tuple powers_i in
                 Array.iteri (fun j_c powers_j ->
-                    let j_c = (Contracted_shell.index basis.(j)) + j_c  in
+                    let j_c = (Contracted_shell.index shell.(j)) + j_c + 1 in
                     let xj = to_int_tuple powers_j in
                     Array.iteri (fun k_c powers_k ->
-                        let k_c = (Contracted_shell.index basis.(k)) + k_c  in 
+                        let k_c = (Contracted_shell.index shell.(k)) + k_c + 1 in 
                         let xk = to_int_tuple powers_k in
                         Array.iteri (fun l_c powers_l ->
-                            let l_c = (Contracted_shell.index basis.(l)) + l_c  in
+                            let l_c = (Contracted_shell.index shell.(l)) + l_c + 1 in
                             let xl = to_int_tuple powers_l in
                             let key = 
                               if swap then
@@ -203,36 +189,36 @@ let to_file ~filename basis =
                               )
                             else
                               out := !out + 1;
-                          ) (Contracted_shell.powers basis.(l)) 
-                      ) (Contracted_shell.powers basis.(k))
-                  ) (Contracted_shell.powers basis.(j))
-              ) (Contracted_shell.powers basis.(i));
-            );
+                          ) (Contracted_shell.powers shell.(l)) 
+                      ) (Contracted_shell.powers shell.(k))
+                  ) (Contracted_shell.powers shell.(j))
+              ) (Contracted_shell.powers shell.(i));
             with NullIntegral ->  ()
           done;
         done;
       with NullIntegral -> ()
-      ) unique_shell_pairs;
-      (*
     done;
   done;
-  *)
-
-  Printf.printf "Computed %d non-zero ERIs in %f seconds\n" !inn (Unix.gettimeofday () -. t0);
-
-  (* Print ERIs *)
-  for i_c=1 to (Genarray.nth_dim eri_array 0) do
-   for j_c=1 to (Genarray.nth_dim eri_array 2) do
-    for k_c=1 to (Genarray.nth_dim eri_array 1) do
-     for l_c=1 to (Genarray.nth_dim eri_array 3) do
-      let value = eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} in
-      if (abs_float value > cutoff) then
-        Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i_c k_c j_c l_c value;
-     done;
-    done;
-   done;
-  done;
   Printf.printf "In: %d  Out:%d\n" !inn !out ;
+  Printf.printf "Computed ERIs in %f seconds\n" (Unix.gettimeofday () -. t0);
+  eri_array
+
+
+(** Write all integrals to a file with the <ij|kl> convention *)
+let to_file ~filename eri_array =
+  let oc = open_out filename in
+  (* Print ERIs *)
+  for l_c=1 to (Genarray.nth_dim eri_array 3) do
+    for k_c=1 to l_c do
+      for j_c=1 to l_c do
+        for i_c=1 to k_c do
+          let value = eri_array.{i_c,j_c,k_c,l_c} in
+          if (abs_float value > cutoff) then
+            Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i_c j_c k_c l_c value;
+        done;
+      done;
+    done;
+  done;
   close_out oc
 
 (*
diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml
index 151d093..d192137 100644
--- a/Basis/KinInt.ml
+++ b/Basis/KinInt.ml
@@ -1,5 +1,9 @@
 open Util
 open Constants
+open Lacaml.D
+
+type t = Mat.t
+
 
 (** Computes all the kinetic integrals of the contracted shell pair *)
 let contracted_class shell_a shell_b : float Zmap.t = 
@@ -94,48 +98,65 @@ let contracted_class shell_a shell_b : float Zmap.t =
   result
 
 
-
-
-
-(** Write all kinetic integrals to a file *)
-let to_file ~filename basis =
-  let to_int_tuple x = 
+(** Create kinetic energy matrix *)
+let of_basis basis =
+  let to_int_tuple x =
     let open Zkey in
     match to_int_tuple Kind_3 x with
     | Three x -> x
     | _ -> assert false
   in
 
-  let oc = open_out filename in
-  for i=0 to (Array.length basis) - 1 do
-    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
-    for j=0 to i do
+  let n =
+    Basis.size basis
+  and shell =
+    Basis.contracted_shells basis
+  in
+
+  let result = Mat.create n n in
+  for j=0 to (Array.length shell) - 1 do
+    for i=0 to j do
       (* Compute all the integrals of the class *)
       let cls =
-        contracted_class basis.(i) basis.(j) 
+        contracted_class shell.(i) shell.(j)
       in
 
-      (* Write the data in the output file *)
-      Array.iteri (fun i_c powers_i ->
-          let i_c = Contracted_shell.index basis.(i) + i_c + 1 in
-          let xi = to_int_tuple powers_i in
-          Array.iteri (fun j_c powers_j ->
-              let j_c = Contracted_shell.index basis.(j) + j_c + 1 in
-              let xj = to_int_tuple powers_j in
-              let key = 
+      Array.iteri (fun j_c powers_j ->
+          let j_c = Contracted_shell.index shell.(j) + j_c + 1 in
+          let xj = to_int_tuple powers_j in
+          Array.iteri (fun i_c powers_i ->
+              let i_c = Contracted_shell.index shell.(i) + i_c + 1 in
+              let xi = to_int_tuple powers_i in
+              let key =
                 Zkey.of_int_tuple (Zkey.Six (xi,xj))
               in
-              let value = 
-                try
-                  Zmap.find cls key 
-                with Not_found -> failwith "Bug in Kinetic integrals"
+              let value =
+                try Zmap.find cls key
+                with Not_found -> failwith "Bug in kinetic integrals"
               in
-              if (abs_float value > cutoff) then
-                Printf.fprintf oc "%4d %4d %20.12e\n"
-                  i_c j_c value
-            ) (Contracted_shell.powers basis.(j))
-        ) (Contracted_shell.powers basis.(i));
+              result.{i_c,j_c} <- value
+            ) (Contracted_shell.powers shell.(i));
+        ) (Contracted_shell.powers shell.(j))
+    done;
+  done;
+  Mat.detri result;
+  result
+
+
+(** Write all kinetic integrals to a file *)
+let to_file ~filename kinetic =
+
+  let oc = open_out filename in
+  let n =
+    Mat.dim1 kinetic
+  in
+
+  for j=1 to n  do
+    for i=1 to j do
+      if (abs_float kinetic.{i,j} > cutoff) then
+        Printf.fprintf oc "%4d %4d %20.12e\n" i j kinetic.{i,j}
     done;
   done;
   close_out oc
 
+
diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml
index d14c934..d28d9bf 100644
--- a/Basis/NucInt.ml
+++ b/Basis/NucInt.ml
@@ -2,7 +2,9 @@
 
 open Util
 open Constants
-open Bigarray
+open Lacaml.D
+
+type t = Mat.t
 
 (** (0|0)^m : Fundamental electron-nucleus repulsion integral
     $ \int \phi_p(r1) 1/r_{C} dr_1 $
@@ -32,7 +34,8 @@ let contracted_class_shell_pair shell_p geometry: float Zmap.t =
 let cutoff2 = cutoff *. cutoff
 exception NullIntegral
 
-let to_file ~filename basis geometry =
+
+let of_basis_nuclei basis nuclei =
   let to_int_tuple x = 
     let open Zkey in
     match to_int_tuple Kind_3 x with
@@ -40,32 +43,22 @@ let to_file ~filename basis geometry =
     | _ -> assert false
   in
 
-  let oc = open_out filename in
+  let n =
+    Basis.size basis
+  and shell =
+    Basis.contracted_shells basis
+  in
+
+  let eni_array = Mat.create n n in
 
   (* Pre-compute all shell pairs *)
   let shell_pairs =
     Array.mapi (fun i shell_a -> Array.map (fun shell_b -> 
-        ContractedShellPair.create shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis
+        ContractedShellPair.create shell_a shell_b) (Array.sub shell 0 (i+1)) ) shell
   in
-  Printf.printf "%d shells\n" (Array.length basis);
-
-  let eni_array = 
-    let n = ref 0 in
-    for i=0 to (Array.length basis) - 1 do
-      n := !n + (Array.length (Contracted_shell.powers basis.(i)))
-    done;
-    let n = !n in
-    Array2.create Float64 c_layout n  n
-  in
-  Array2.fill eni_array 0.;
 
   (* Compute Integrals *)
-  let t0 = Unix.gettimeofday () in
-
-  let inn = ref 0 and out = ref 0 in
-
-  for i=0 to (Array.length basis) - 1 do
-    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
+  for i=0 to (Array.length shell) - 1 do
     for j=0 to i do
       let
         shell_p = shell_pairs.(i).(j)
@@ -73,15 +66,15 @@ let to_file ~filename basis geometry =
 
       (* Compute all the integrals of the class *)
       let cls =
-        contracted_class_shell_pair shell_p geometry
+        contracted_class_shell_pair shell_p nuclei
       in
 
       (* Write the data in the output file *)
       Array.iteri (fun i_c powers_i ->
-          let i_c = (Contracted_shell.index basis.(i)) + i_c + 1 in
+          let i_c = (Contracted_shell.index shell.(i)) + i_c + 1 in
           let xi = to_int_tuple powers_i in
           Array.iteri (fun j_c powers_j ->
-              let j_c = (Contracted_shell.index basis.(j)) + j_c + 1 in
+              let j_c = (Contracted_shell.index shell.(j)) + j_c + 1 in
               let xj = to_int_tuple powers_j in
               let key = 
                 Zkey.of_int_tuple (Zkey.Six (xi,xj))
@@ -89,27 +82,25 @@ let to_file ~filename basis geometry =
               let value = 
                 Zmap.find cls key 
               in
-              if (abs_float value > cutoff) then
-                (inn := !inn + 1;
-                 eni_array.{(i_c-1),(j_c-1)} <- value;
-                )
-              else
-                out := !out + 1;
-            ) (Contracted_shell.powers basis.(j))
-        ) (Contracted_shell.powers basis.(i));
+              eni_array.{j_c,i_c} <- value;
+            ) (Contracted_shell.powers shell.(j))
+        ) (Contracted_shell.powers shell.(i));
     done;
   done;
+  Mat.detri eni_array;
+  eni_array
 
-  Printf.printf "Computed %d non-zero ENIs in %f seconds\n" !inn (Unix.gettimeofday () -. t0);
 
-  (* Print ENIs *)
-  for i_c=1 to (Array2.dim1 eni_array) do
-    for j_c=1 to (Array2.dim2 eni_array) do
-      let value = eni_array.{(i_c-1),(j_c-1)} in
+let to_file ~filename eni_array =
+  let n = Mat.dim1 eni_array in
+  let oc = open_out filename in
+
+  for j=1 to n do
+    for i=1 to j do
+      let value = eni_array.{i,j} in
       if (value <> 0.) then
-        Printf.fprintf oc " %5d %5d %20.15f\n" i_c j_c value;
+        Printf.fprintf oc " %5d %5d %20.15f\n" i j value;
     done;
   done;
-  Printf.printf "In: %d  Out:%d\n" !inn !out ;
   close_out oc
 
diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml
index 0e73d50..740d40c 100644
--- a/Basis/Overlap.ml
+++ b/Basis/Overlap.ml
@@ -1,5 +1,9 @@
 open Util
 open Constants
+open Lacaml.D
+
+type t = Mat.t
+
 
 (** Computes all the overlap integrals of the contracted shell pair *)
 let contracted_class shell_a shell_b : float Zmap.t = 
@@ -66,11 +70,8 @@ let contracted_class shell_a shell_b : float Zmap.t =
   result
 
 
-
-
-
-(** Write all overlap integrals to a file *)
-let to_file ~filename basis =
+(** Create overlap matrix *)
+let of_basis basis =
   let to_int_tuple x = 
     let open Zkey in
     match to_int_tuple Kind_3 x with
@@ -78,35 +79,56 @@ let to_file ~filename basis =
     | _ -> assert false
   in
 
-  let oc = open_out filename in
-  for i=0 to (Array.length basis) - 1 do
-    print_int (Contracted_shell.index basis.(i)) ; print_newline ();
-    for j=0 to i do
+  let n =
+    Basis.size basis
+  and shell =
+    Basis.contracted_shells basis
+  in
+
+  let result = Mat.create n n in
+  for j=0 to (Array.length shell) - 1 do
+    for i=0 to j do
       (* Compute all the integrals of the class *)
       let cls =
-        contracted_class basis.(i) basis.(j) 
+        contracted_class shell.(i) shell.(j) 
       in
 
-      (* Write the data in the output file *)
-      Array.iteri (fun i_c powers_i ->
-          let i_c = Contracted_shell.index basis.(i) + i_c + 1 in
-          let xi = to_int_tuple powers_i in
-          Array.iteri (fun j_c powers_j ->
-              let j_c = Contracted_shell.index basis.(j) + j_c + 1 in
-              let xj = to_int_tuple powers_j in
+      Array.iteri (fun j_c powers_j ->
+          let j_c = Contracted_shell.index shell.(j) + j_c + 1 in
+          let xj = to_int_tuple powers_j in
+          Array.iteri (fun i_c powers_i ->
+              let i_c = Contracted_shell.index shell.(i) + i_c + 1 in
+              let xi = to_int_tuple powers_i in
               let key = 
                 Zkey.of_int_tuple (Zkey.Six (xi,xj))
               in
               let value = 
-                try
-                  Zmap.find cls key 
+                try Zmap.find cls key 
                 with Not_found -> failwith "Bug in overlap integrals"
               in
-              if (abs_float value > cutoff) then
-                Printf.fprintf oc "%4d %4d %20.12e\n"
-                  i_c j_c value
-            ) (Contracted_shell.powers basis.(j))
-        ) (Contracted_shell.powers basis.(i));
+              result.{i_c,j_c} <- value
+            ) (Contracted_shell.powers shell.(i));
+        ) (Contracted_shell.powers shell.(j))
+    done;
+  done;
+  Mat.detri result;
+  result
+
+
+
+
+(** Write all overlap integrals to a file *)
+let to_file ~filename overlap =
+
+  let oc = open_out filename in
+  let n =
+    Mat.dim1 overlap 
+  in
+
+  for j=1 to n  do
+    for i=1 to j do
+      if (abs_float overlap.{i,j} > cutoff) then
+        Printf.fprintf oc "%4d %4d %20.12e\n" i j overlap.{i,j}
     done;
   done;
   close_out oc
diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml
index e9d5712..c1babb6 100644
--- a/Nuclei/Nuclei.ml
+++ b/Nuclei/Nuclei.ml
@@ -23,6 +23,7 @@ let of_zmt_file ~filename =
   |> Zmatrix.to_xyz
   |> Array.map (fun (e,x,y,z) -> (e, Coordinate.of_3_floats x y z ))
     
+
 let to_string atoms =
     "
                 Nuclear Coordinates (Angstrom)
@@ -49,14 +50,7 @@ let to_string atoms =
 "
 
 
-let file : string option ref = ref None
-
-let set_file f =
-  file := Some f
-
-let nuclei = lazy(
-  match !file with
-  | None -> failwith "coordinate file not defined"
-  | Some filename -> of_xyz_file ~filename
-)   
+let of_filename ~filename =
+    of_xyz_file filename
+  
 
diff --git a/README.rst b/README.rst
index 233261b..bd8b387 100644
--- a/README.rst
+++ b/README.rst
@@ -2,8 +2,8 @@ Requirements
 ------------
 
 * gmp : GNU Multiple Precision arithmetic library
-* zarith : Arbitrary-precision integers
 * BLAS/LAPACK : Linear algebra
+* Zarith : Arbitrary-precision integers
 * LaCaml : LAPACK OCaml interface
 * SklMl : Parallel skeletons for OCaml
 
diff --git a/_tags b/_tags
index 870e647..e3a59c1 100644
--- a/_tags
+++ b/_tags
@@ -1,3 +1,3 @@
-true: package(str,unix,bigarray,zarith)
+true: package(str,unix,bigarray,zarith,lacaml)
 <*.byte>  : linkdep(Utils/math_functions.o), custom
 <*.native>: linkdep(Utils/math_functions.o)
diff --git a/run_integrals.ml b/run_integrals.ml
index f7bdc6e..c165659 100644
--- a/run_integrals.ml
+++ b/run_integrals.ml
@@ -1,10 +1,12 @@
-let out_file   : string option ref = ref None
+let out_file    : string option ref = ref None
+let basis_file  : string option ref = ref None
+let nuclei_file : string option ref = ref None
 
 
 let speclist = [
-  ( "-b" , Arg.String Basis.set_file ,
+  ( "-b" , Arg.String (fun x -> basis_file := Some x),
     "File containing the atomic basis set") ;
-  ( "-c" , Arg.String Nuclei.set_file ,
+  ( "-c" , Arg.String (fun x -> nuclei_file := Some x),
     "File containing the nuclear coordinates") ;
   ( "-o" , Arg.String (fun x -> out_file := Some x) ,
     "Output file") ;
@@ -19,18 +21,31 @@ let run ~out =
     match out with 
     | None -> raise (Invalid_argument "Output file should be specified with -o")
     | Some x -> x
+  and basis_file = 
+    match !basis_file with
+    | None -> raise (Invalid_argument "Basis set file should be specified with -b")
+    | Some x -> x
+  and nuclei_file = 
+    match !nuclei_file with
+    | None -> raise (Invalid_argument "Basis set file should be specified with -c")
+    | Some x -> x
   in
 
-  let nuclei = Lazy.force Nuclei.nuclei in
-  print_endline @@ Nuclei.to_string nuclei;
+  let s =
+    Simulation.of_filenames ~nuclei:nuclei_file ~basis:basis_file 
+  in
 
-  let basis = Lazy.force Basis.basis in
-  print_endline @@ Basis.to_string basis;
+  print_endline @@ Nuclei.to_string s.Simulation.nuclei;
+  print_endline @@ Basis.to_string  s.Simulation.basis;
 
-  Overlap.to_file ~filename:(out_file^".overlap") basis;
-  NucInt.to_file ~filename:(out_file^".nuc") basis nuclei;
-  KinInt.to_file ~filename:(out_file^".kin") basis;
-  ERI.to_file ~filename:(out_file^".eri") basis
+  let overlap  = Lazy.force s.Simulation.overlap  in
+  let eN_ints  = Lazy.force s.Simulation.eN_ints  in
+  let kin_ints = Lazy.force s.Simulation.kin_ints in
+  let ee_ints  = Lazy.force s.Simulation.ee_ints  in
+  Overlap.to_file ~filename:(out_file^".overlap") overlap;
+  NucInt.to_file ~filename:(out_file^".nuc") eN_ints;
+  KinInt.to_file ~filename:(out_file^".kin") kin_ints;
+  ERI.to_file ~filename:(out_file^".eri") ee_ints
 
 
 let () =