diff --git a/.gitignore b/.gitignore
index 81626aa..3ac40d5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@ _build
 .merlin
 *.install
 qcaml.opam
+#*#
diff --git a/ao/lib/basis.ml b/ao/lib/basis.ml
index 808aeda..7ded4c8 100644
--- a/ao/lib/basis.ml
+++ b/ao/lib/basis.ml
@@ -1,11 +1,8 @@
 open Linear_algebra
-    
-type basis =
-  | Unknown
-  | Gaussian of Basis_gaussian.t
+open Common
 
 type t =
-  { ao_basis : basis  ;
+  { ao_basis : Basis_poly.t  ;
     cartesian : bool
   }
 
@@ -19,25 +16,24 @@ let of_nuclei_and_basis_filename ?(kind=`Gaussian) ?operators ?(cartesian=false)
         Gaussian.Basis.of_nuclei_and_basis_filename  ~nuclei filename
       in
       let ao_basis = 
-        Gaussian (Basis_gaussian.make ~basis ?operators ~cartesian nuclei )
+        Basis_poly.Gaussian (Basis_gaussian.make ~basis ?operators ~cartesian nuclei )
       in
       { ao_basis ; cartesian }
   | _ -> failwith "of_nuclei_and_basis_filename needs to be called with `Gaussian"
 
-let not_implemented () =
-  failwith "Not implemented"
+let not_implemented = Util.not_implemented 
   
 let ao_basis t = t.ao_basis
                     
 let size t = 
   match t.ao_basis with
-  | Gaussian b -> Basis_gaussian.size b
+  | Basis_poly.Gaussian b -> Basis_gaussian.size b
   | _ -> not_implemented ()
 
 let overlap t =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.overlap b
+    | Basis_poly.Gaussian b -> Basis_gaussian.overlap b
     | _ -> not_implemented ()
   end
   |> Matrix.relabel
@@ -45,15 +41,18 @@ let overlap t =
 let multipole t =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.multipole b
+    | Basis_poly.Gaussian b ->
+        let m = Basis_gaussian.multipole b in
+        fun s ->
+          Gaussian_integrals.Multipole.matrix m s
+          |> Matrix.relabel
     | _ -> not_implemented ()
   end
-  |> Array.map Matrix.relabel
        
 let ortho t =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.ortho b
+    | Basis_poly.Gaussian b -> Basis_gaussian.ortho b
     | _ -> not_implemented ()
   end
   |> Matrix.relabel
@@ -61,7 +60,7 @@ let ortho t =
 let eN_ints t =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.eN_ints b
+    | Basis_poly.Gaussian b -> Basis_gaussian.eN_ints b
     | _ -> not_implemented ()
   end
   |> Matrix.relabel
@@ -69,7 +68,7 @@ let eN_ints t =
 let kin_ints t =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.kin_ints b
+    | Basis_poly.Gaussian b -> Basis_gaussian.kin_ints b
     | _ -> not_implemented ()
   end
   |> Matrix.relabel
@@ -77,7 +76,7 @@ let kin_ints t =
 let ee_ints t =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.ee_ints b
+    | Basis_poly.Gaussian b -> Basis_gaussian.ee_ints b
     | _ -> not_implemented ()
   end
   |> Four_idx_storage.relabel
@@ -85,7 +84,7 @@ let ee_ints t =
 let ee_lr_ints t =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.ee_lr_ints b
+    | Basis_poly.Gaussian b -> Basis_gaussian.ee_lr_ints b
     | _ -> not_implemented ()
   end
   |> Four_idx_storage.relabel
@@ -93,7 +92,7 @@ let ee_lr_ints t =
 let f12_ints t =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.f12_ints b
+    | Basis_poly.Gaussian b -> Basis_gaussian.f12_ints b
     | _ -> not_implemented ()
   end
   |> Four_idx_storage.relabel
@@ -101,7 +100,7 @@ let f12_ints t =
 let f12_over_r12_ints t =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.f12_over_r12_ints b
+    | Basis_poly.Gaussian b -> Basis_gaussian.f12_over_r12_ints b
     | _ -> not_implemented ()
   end
   |> Four_idx_storage.relabel
@@ -112,7 +111,7 @@ let cartesian t = t.cartesian
 let values t point =
   begin
     match t.ao_basis with
-    | Gaussian b -> Basis_gaussian.values b point
+    | Basis_poly.Gaussian b -> Basis_gaussian.values b point
     | _ -> not_implemented ()
   end
   |> Vector.relabel
diff --git a/ao/lib/basis.mli b/ao/lib/basis.mli
index 38021b0..c4e38ce 100644
--- a/ao/lib/basis.mli
+++ b/ao/lib/basis.mli
@@ -5,10 +5,6 @@ open Particles
 open Operators
 open Linear_algebra
 
-type basis =
-  | Unknown
-  | Gaussian of Basis_gaussian.t
-
 type t 
 type ao = Ao_dim.t
 
@@ -17,13 +13,13 @@ type ao = Ao_dim.t
 val size : t -> int
 (** Number of atomic orbitals in the AO basis set *)
   
-val ao_basis : t -> basis
+val ao_basis : t -> Basis_poly.t
 (** One-electron basis set *)
 
 val overlap : t -> (ao,ao) Matrix.t
 (** Overlap matrix *)
 
-val multipole : t -> (ao,ao) Matrix.t array
+val multipole : t -> string -> (ao,ao) Matrix.t
 (** Multipole matrices *)
 
 val ortho : t -> (ao,'a) Matrix.t
diff --git a/ao/lib/basis_poly.ml b/ao/lib/basis_poly.ml
new file mode 100644
index 0000000..d393bd1
--- /dev/null
+++ b/ao/lib/basis_poly.ml
@@ -0,0 +1,3 @@
+type t =
+  | Unknown
+  | Gaussian of Basis_gaussian.t
diff --git a/ao/lib/basis_poly.mli b/ao/lib/basis_poly.mli
new file mode 100644
index 0000000..10a5770
--- /dev/null
+++ b/ao/lib/basis_poly.mli
@@ -0,0 +1,6 @@
+(** Polymorphic type for the basis (Gaussian, Slater, etc) *)
+
+type t =
+  | Unknown
+  | Gaussian of Basis_gaussian.t
+    
diff --git a/common/lib/util.ml b/common/lib/util.ml
index a225d1b..103501b 100644
--- a/common/lib/util.ml
+++ b/common/lib/util.ml
@@ -43,6 +43,8 @@ let () =
   Sys.set_signal Sys.sigint (Sys.Signal_handle f) 
   ;;
 
+let not_implemented () = 
+   failwith "Not implemented"
 
 let memo_float_of_int =
   Array.init 64 float_of_int
diff --git a/common/lib/util.mli b/common/lib/util.mli
index 316bb64..fb02c2b 100644
--- a/common/lib/util.mli
+++ b/common/lib/util.mli
@@ -41,6 +41,9 @@ val leadz : int64 -> int
 
 (** {2 General functions} *)
 
+val not_implemented : unit -> 'a
+(** Fails with error message if some functionality is not implemented. *)
+  
 val fact : int -> float
 (** Factorial function.
   @raise Invalid_argument for negative arguments or arguments >100.
diff --git a/examples/dune b/examples/dune
new file mode 100644
index 0000000..accdf4d
--- /dev/null
+++ b/examples/dune
@@ -0,0 +1,6 @@
+(executable
+ (name ex_integrals)
+ (libraries
+  qcaml
+))
+
diff --git a/examples/ex_integrals.ml b/examples/ex_integrals.ml
new file mode 100644
index 0000000..f27405b
--- /dev/null
+++ b/examples/ex_integrals.ml
@@ -0,0 +1,93 @@
+(* [[file:~/QCaml/examples/ex_integrals.org::*Header][Header:1]] *)
+module Command_line = Qcaml.Common.Command_line
+module Util = Qcaml.Common.Util
+open Qcaml.Linear_algebra
+
+let () =
+(* Header:1 ends here *)
+
+(* [[file:~/QCaml/examples/ex_integrals.org::*Definition][Definition:1]] *)
+let open Command_line in
+begin
+  set_header_doc (Sys.argv.(0) ^ " - QuAcK command");
+  set_description_doc "Computes the one- and two-electron integrals on the Gaussian atomic basis set.";
+  set_specs
+    [ { short='b' ; long="basis" ; opt=Mandatory;
+	arg=With_arg "<string>";
+	doc="Name of the file containing the basis set"; } ;
+
+      { short='x' ; long="xyz" ; opt=Mandatory;
+	arg=With_arg "<string>";
+	doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
+
+      { short='u' ; long="range-separation" ; opt=Optional;
+	arg=With_arg "<float>";
+	doc="Range-separation parameter."; } ;
+    ]
+end;
+(* Definition:1 ends here *)
+
+(* [[file:~/QCaml/examples/ex_integrals.org::*Interpretation][Interpretation:1]] *)
+let basis_file  = Util.of_some @@ Command_line.get "basis" in
+let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
+
+let range_separation =
+  match Command_line.get "range-separation" with
+  | None -> None
+  | Some mu -> Some (float_of_string mu)
+in
+
+let operators =
+  match range_separation with
+  | None -> []
+  | Some mu -> [ Qcaml.Operators.Operator.of_range_separation mu ]
+in
+(* Interpretation:1 ends here *)
+
+(* [[file:~/QCaml/examples/ex_integrals.org::*Computation][Computation:1]] *)
+let nuclei =
+  Qcaml.Particles.Nuclei.of_xyz_file nuclei_file
+in
+(* Computation:1 ends here *)
+
+(* [[file:~/QCaml/examples/ex_integrals.org::*Computation][Computation:2]] *)
+let ao_basis =
+  Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian
+    ~operators ~cartesian:true ~nuclei basis_file
+in
+(* Computation:2 ends here *)
+
+(* [[file:~/QCaml/examples/ex_integrals.org::*Computation][Computation:3]] *)
+let overlap   = Qcaml.Ao.Basis.overlap   ao_basis in
+let eN_ints   = Qcaml.Ao.Basis.eN_ints   ao_basis in
+let kin_ints  = Qcaml.Ao.Basis.kin_ints  ao_basis in
+let multipole = Qcaml.Ao.Basis.multipole ao_basis in
+let x_mat     = multipole "x" in
+let y_mat     = multipole "y" in
+let z_mat     = multipole "z" in
+(* Computation:3 ends here *)
+
+(* [[file:~/QCaml/examples/ex_integrals.org::*Computation][Computation:4]] *)
+let ee_ints   = Qcaml.Ao.Basis.ee_ints   ao_basis in
+let lr_ints   = 
+  match range_separation with
+  | Some _mu -> Some (Qcaml.Ao.Basis.ee_lr_ints ao_basis)
+  | None -> None
+in
+(* Computation:4 ends here *)
+
+(* [[file:~/QCaml/examples/ex_integrals.org::*Output][Output:1]] *)
+Matrix.to_file ~filename:"overlap.dat" ~sym:true overlap;
+Matrix.to_file ~filename:"eN.dat" ~sym:true eN_ints;
+Matrix.to_file ~filename:"kinetic.dat" ~sym:true kin_ints;
+Matrix.to_file ~filename:"x.dat" ~sym:true x_mat;
+Matrix.to_file ~filename:"y.dat" ~sym:true y_mat;
+Matrix.to_file ~filename:"z.dat" ~sym:true z_mat;
+(* Output:1 ends here *)
+
+(* [[file:~/QCaml/examples/ex_integrals.org::*Output][Output:2]] *)
+Four_idx_storage.to_file ~filename:"eri.dat" ee_ints;
+match lr_ints with
+| Some integrals -> Four_idx_storage.to_file ~filename:"eri_lr.dat" integrals;
+| None -> ()
+(* Output:2 ends here *)
diff --git a/examples/ex_integrals.org b/examples/ex_integrals.org
new file mode 100644
index 0000000..0ce91cf
--- /dev/null
+++ b/examples/ex_integrals.org
@@ -0,0 +1,130 @@
+#+TITLE: Integrals
+
+#+PROPERTY
+
+In this example, we write a program that reads the geometry of a
+molecule in in =xyz= format and a Gaussian atomic basis set in GAMESS
+format. The output is a set of files containing the one- and two-
+electron integrals.
+
+* Header 
+
+#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+module Command_line = Qcaml.Common.Command_line
+module Util = Qcaml.Common.Util
+open Qcaml.Linear_algebra
+
+let () =
+#+END_SRC
+
+* Command-line arguments
+
+  We use the =Command_line= module to define the following possible
+  arguments:
+  - =-b --basis= : The name of the file containing the basis set
+  - =-x --xyz=   : The name of the file containing the atomic coordinates
+  - =-u --range-separation= : The value of $\mu$, the range-separation
+    parameter in range-separated DFT. If this option is not present,
+    no output file will be generated for the range-separated integrals.
+
+** Definition
+
+   #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let open Command_line in
+begin
+  set_header_doc (Sys.argv.(0) ^ " - QuAcK command");
+  set_description_doc "Computes the one- and two-electron integrals on the Gaussian atomic basis set.";
+  set_specs
+    [ { short='b' ; long="basis" ; opt=Mandatory;
+	arg=With_arg "<string>";
+	doc="Name of the file containing the basis set"; } ;
+
+      { short='x' ; long="xyz" ; opt=Mandatory;
+	arg=With_arg "<string>";
+	doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
+
+      { short='u' ; long="range-separation" ; opt=Optional;
+	arg=With_arg "<float>";
+	doc="Range-separation parameter."; } ;
+    ]
+end;
+   #+END_SRC
+
+** Interpretation
+
+   #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let basis_file  = Util.of_some @@ Command_line.get "basis" in
+let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
+
+let range_separation =
+  match Command_line.get "range-separation" with
+  | None -> None
+  | Some mu -> Some (float_of_string mu)
+in
+
+let operators =
+  match range_separation with
+  | None -> []
+  | Some mu -> [ Qcaml.Operators.Operator.of_range_separation mu ]
+in
+   #+END_SRC
+
+* Computation
+  We first read the =xyz= file to create a molecule:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let nuclei =
+  Qcaml.Particles.Nuclei.of_xyz_file nuclei_file
+in
+  #+END_SRC
+
+  Then we create an Gaussian AO basis using the atomic coordinates,
+  and we optionally introduce the range-separation parameter via the
+  =operators=:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let ao_basis =
+  Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian
+    ~operators ~cartesian:true ~nuclei basis_file
+in
+  #+END_SRC
+
+  We compute the required one-electron integrals:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let overlap   = Qcaml.Ao.Basis.overlap   ao_basis in
+let eN_ints   = Qcaml.Ao.Basis.eN_ints   ao_basis in
+let kin_ints  = Qcaml.Ao.Basis.kin_ints  ao_basis in
+let multipole = Qcaml.Ao.Basis.multipole ao_basis in
+let x_mat     = multipole "x" in
+let y_mat     = multipole "y" in
+let z_mat     = multipole "z" in
+  #+END_SRC
+
+  and the two-electron integrals (1/r and long range):
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let ee_ints   = Qcaml.Ao.Basis.ee_ints   ao_basis in
+let lr_ints   = 
+  match range_separation with
+  | Some _mu -> Some (Qcaml.Ao.Basis.ee_lr_ints ao_basis)
+  | None -> None
+in
+  #+END_SRC
+
+* Output
+
+  We write the one-electron integrals:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+Matrix.to_file ~filename:"overlap.dat" ~sym:true overlap;
+Matrix.to_file ~filename:"eN.dat" ~sym:true eN_ints;
+Matrix.to_file ~filename:"kinetic.dat" ~sym:true kin_ints;
+Matrix.to_file ~filename:"x.dat" ~sym:true x_mat;
+Matrix.to_file ~filename:"y.dat" ~sym:true y_mat;
+Matrix.to_file ~filename:"z.dat" ~sym:true z_mat;
+  #+END_SRC
+
+  and the the two-electron integrals:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+Four_idx_storage.to_file ~filename:"eri.dat" ee_ints;
+match lr_ints with
+| Some integrals -> Four_idx_storage.to_file ~filename:"eri_lr.dat" integrals;
+| None -> ()
+  #+END_SRC
+  
diff --git a/gaussian_integrals/lib/electron_nucleus.ml b/gaussian_integrals/lib/electron_nucleus.ml
index 74eaaea..db7f027 100644
--- a/gaussian_integrals/lib/electron_nucleus.ml
+++ b/gaussian_integrals/lib/electron_nucleus.ml
@@ -91,21 +91,6 @@ let of_basis_nuclei ~basis nuclei =
   Matrix.detri_inplace eni_array;
   eni_array
 
-
-let to_file ~filename eni_array =
-  let n = Matrix.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 j value;
-    done;
-  done;
-  close_out oc
-
-
 let of_basis _basis =
   invalid_arg "of_basis_nuclei should be called for NucInt"
 
diff --git a/gaussian_integrals/lib/kinetic.ml b/gaussian_integrals/lib/kinetic.ml
index 6b045da..214fd6e 100644
--- a/gaussian_integrals/lib/kinetic.ml
+++ b/gaussian_integrals/lib/kinetic.ml
@@ -166,20 +166,3 @@ let of_basis basis =
 let of_basis_pair _first_basis _second_basis =
   failwith "Not implemented"
 
-(** Write all kinetic integrals to a file *)
-let to_file ~filename kinetic =
-
-  let oc = open_out filename in
-  let n =
-    Matrix.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/gaussian_integrals/lib/matrix_on_basis.mli b/gaussian_integrals/lib/matrix_on_basis.mli
index 2ad40bd..6288c4e 100644
--- a/gaussian_integrals/lib/matrix_on_basis.mli
+++ b/gaussian_integrals/lib/matrix_on_basis.mli
@@ -10,6 +10,3 @@ val of_basis : Basis.t -> t
 
 val of_basis_pair : Basis.t -> Basis.t -> t
 (** Build from a pair of {!Basis.t}. *)
-
-val to_file : filename:string -> t -> unit
-(** Write the integrals in a file. *)
diff --git a/gaussian_integrals/lib/multipole.ml b/gaussian_integrals/lib/multipole.ml
index 83cfe40..ba8a66e 100644
--- a/gaussian_integrals/lib/multipole.ml
+++ b/gaussian_integrals/lib/multipole.ml
@@ -16,22 +16,23 @@ module Csp = Contracted_shell_pair
 module Po  = Powers
 module Psp = Primitive_shell_pair
 
-let matrix_x  t =  t.(0)
-let matrix_y  t =  t.(1)
-let matrix_z  t =  t.(2)
-let matrix_x2 t =  t.(3)
-let matrix_y2 t =  t.(4)
-let matrix_z2 t =  t.(5)
-let matrix_xy t =  t.(6)
-let matrix_xz t =  t.(8)
-let matrix_yz t =  t.(7)
-let matrix_x3 t =  t.(9)
-let matrix_y3 t =  t.(10)
-let matrix_z3 t =  t.(11)
-let matrix_x4 t =  t.(12)
-let matrix_y4 t =  t.(13)
-let matrix_z4 t =  t.(14)
-
+let matrix t = function
+  | "x" -> t.(0)
+  | "y" -> t.(1)
+  | "z" -> t.(2)
+  | "x2" -> t.(3)
+  | "y2" -> t.(4)
+  | "z2" -> t.(5)
+  | "xy" -> t.(6)
+  | "xz" -> t.(8)
+  | "yz" -> t.(7)
+  | "x3" -> t.(9)
+  | "y3" -> t.(10)
+  | "z3" -> t.(11)
+  | "x4" -> t.(12)
+  | "y4" -> t.(13)
+  | "z4" -> t.(14)
+  | _ -> Util.not_implemented ()
 
 
 
diff --git a/gaussian_integrals/lib/multipole.mli b/gaussian_integrals/lib/multipole.mli
index f897ccc..69dea0d 100644
--- a/gaussian_integrals/lib/multipole.mli
+++ b/gaussian_integrals/lib/multipole.mli
@@ -14,50 +14,10 @@ open Gaussian
 
 type t = (Basis.t, Basis.t) Matrix.t array
 
-val matrix_x : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | x | \chi_j \rangle $$ %} *)
-
-val matrix_y : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | y | \chi_j \rangle $$ %} *)
-
-val matrix_z : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | z | \chi_j \rangle $$ %} *)
-
-val matrix_x2 : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %} *)
-
-val matrix_xy : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | xy | \chi_j \rangle $$ %} *)
-
-val matrix_yz : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | yz | \chi_j \rangle $$ %} *)
-
-val matrix_xz : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | xz | \chi_j \rangle $$ %} *)
-
-val matrix_y2 : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %} *)
-
-val matrix_z2 : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *)
-
-val matrix_x3 : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | x^3 | \chi_j \rangle $$ %} *)
-
-val matrix_y3 : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | y^3 | \chi_j \rangle $$ %} *)
-
-val matrix_z3 : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | z^3 | \chi_j \rangle $$ %} *)
-
-val matrix_x4 : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | x^4 | \chi_j \rangle $$ %} *)
-
-val matrix_y4 : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | y^4 | \chi_j \rangle $$ %} *)
-
-val matrix_z4 : t -> (Basis.t, Basis.t) Matrix.t
-(** {% $$ \langle \chi_i | z^4 | \chi_j \rangle $$ %} *)
+val matrix : t -> string -> (Basis.t, Basis.t) Matrix.t
+(** Returns the requested matrix. Choices:
+[x, y, z, x2, xy, yz, xz, y2, z2, z3, y3, z3, x4, y4, z4]
+*)
 
 val of_basis : Basis.t -> t
 
diff --git a/gaussian_integrals/lib/overlap.ml b/gaussian_integrals/lib/overlap.ml
index bd7f687..6f13fae 100644
--- a/gaussian_integrals/lib/overlap.ml
+++ b/gaussian_integrals/lib/overlap.ml
@@ -174,21 +174,3 @@ let of_basis_pair first_basis second_basis =
 
 
 
-(** Write all overlap integrals to a file *)
-let to_file ~filename overlap =
-
-  let oc = open_out filename in
-  let n =
-    Matrix.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/linear_algebra/lib/matrix.ml b/linear_algebra/lib/matrix.ml
index 8a4927e..f01f700 100644
--- a/linear_algebra/lib/matrix.ml
+++ b/linear_algebra/lib/matrix.ml
@@ -372,8 +372,37 @@ let qr a =
   let q = result in
   q, r
 
+  
 let (%:) t (i,j) = t.{i,j}
 
+
+let to_file ~filename ?(sym=false) ?(cutoff=0.) t =
+
+  let oc = open_out filename in
+  let n = Mat.dim1 t in
+  let m = Mat.dim2 t in
+
+  if sym then 
+    for j=1 to n  do
+      for i=1 to j do
+        if (abs_float (t.{i,j}) > cutoff) then
+          Printf.fprintf oc "%4d %4d %20.12e\n" i j (t.{i,j})
+      done;
+    done
+  else
+    for j=1 to n  do
+      for i=1 to m do
+        if (abs_float (t.{i,j}) > cutoff) then
+          Printf.fprintf oc "%4d %4d %20.12e\n" i j (t.{i,j})
+      done;
+    done;
+  close_out oc
+
+
+
+  
+  
+
 let set  t i j v = t.{i,j} <- v
 
 
diff --git a/linear_algebra/lib/matrix.mli b/linear_algebra/lib/matrix.mli
index 9b3d92e..6d4011e 100644
--- a/linear_algebra/lib/matrix.mli
+++ b/linear_algebra/lib/matrix.mli
@@ -319,5 +319,8 @@ val (%:) : ('a,'b) t -> (int*int) -> float
 val set : ('a,'b) t -> int -> int -> float -> unit
 (** [set t i j v] sets the (i,j)-th element to v *)
 
+val to_file : filename:string -> ?sym:bool -> ?cutoff:float -> ('a,'b) t -> unit
+(** Write the matrix to a file. *)
+  
 val pp : Format.formatter -> ('a,'b) t -> unit