diff --git a/examples/dune b/examples/dune
index accdf4d..3e5d340 100644
--- a/examples/dune
+++ b/examples/dune
@@ -1,6 +1,7 @@
-(executable
- (name ex_integrals)
- (libraries
-  qcaml
-))
+(executables
+  (names
+        ex_integrals
+        ex_hartree_fock
+)
+(libraries qcaml))
 
diff --git a/examples/ex_hartree_fock.ml b/examples/ex_hartree_fock.ml
new file mode 100644
index 0000000..76fc8cf
--- /dev/null
+++ b/examples/ex_hartree_fock.ml
@@ -0,0 +1,76 @@
+(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Header][Header:1]] *)
+module Command_line = Qcaml.Common.Command_line
+module Util = Qcaml.Common.Util
+
+let () =
+(* Header:1 ends here *)
+
+(* [[file:~/QCaml/examples/ex_hartree_fock.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 hartree_fock 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='m' ; long="multiplicity" ; opt=Optional;
+        arg=With_arg "<int>";
+        doc="Spin multiplicity (2S+1). Default is singlet"; } ;
+      
+      { short='c' ; long="charge" ; opt=Optional;
+        arg=With_arg "<int>";
+        doc="Total charge of the molecule. Specify negative charges with 'm' instead of the minus sign, for example m1 instead of -1. Default is 0"; } ;
+      
+    ]
+end;
+(* Definition:1 ends here *)
+
+(* [[file:~/QCaml/examples/ex_hartree_fock.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 charge =
+  match Command_line.get "charge" with
+  | Some x -> ( if x.[0] = 'm' then
+                  ~- (int_of_string (String.sub x 1 (String.length x - 1)))
+                else
+                  int_of_string x )
+  | None   -> 0
+in
+    
+let multiplicity =
+  match Command_line.get "multiplicity" with
+  | Some x -> int_of_string x
+  | None -> 1
+in
+(* Interpretation:1 ends here *)
+
+(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:1]] *)
+let nuclei =
+  Qcaml.Particles.Nuclei.of_xyz_file nuclei_file
+in
+(* Computation:1 ends here *)
+
+(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:2]] *)
+let ao_basis =
+  Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file
+in
+(* Computation:2 ends here *)
+
+(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:3]] *)
+let simulation = Qcaml.Simulation.make ~charge ~multiplicity ~nuclei ao_basis in
+(* Computation:3 ends here *)
+
+(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:4]] *)
+let hf = Qcaml.Mo.Hartree_fock.make ~guess:`Huckel simulation in
+(* Computation:4 ends here *)
+
+(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Output][Output:1]] *)
+Format.printf "@[%a@]" (Mo.Hartree_fock.pp) hf
+(* Output:1 ends here *)
diff --git a/examples/ex_hartree_fock.org b/examples/ex_hartree_fock.org
new file mode 100644
index 0000000..3887973
--- /dev/null
+++ b/examples/ex_hartree_fock.org
@@ -0,0 +1,109 @@
+#+TITLE: Hartree-Fock
+
+#+PROPERTY
+
+In this example, we write a program that makes a Hartree-Fock
+caculation.  The molecule is read in =xyz= format and a Gaussian
+atomic basis set in GAMESS format. 
+
+* Header 
+
+#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml  
+module Command_line = Qcaml.Common.Command_line
+module Util = Qcaml.Common.Util
+
+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
+  - =-c --charge=  : The charge of the molecule
+  - =-m --multiplicity= : The spin multiplicity
+
+** Definition
+
+   #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml  
+let open Command_line in
+begin
+  set_header_doc (Sys.argv.(0) ^ " - QuAcK command");
+  set_description_doc "Computes the one- and two-electron hartree_fock 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='m' ; long="multiplicity" ; opt=Optional;
+        arg=With_arg "<int>";
+        doc="Spin multiplicity (2S+1). Default is singlet"; } ;
+      
+      { short='c' ; long="charge" ; opt=Optional;
+        arg=With_arg "<int>";
+        doc="Total charge of the molecule. Specify negative charges with 'm' instead of the minus sign, for example m1 instead of -1. Default is 0"; } ;
+      
+    ]
+end;
+   #+END_SRC
+
+** Interpretation
+
+   #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml  
+let basis_file  = Util.of_some @@ Command_line.get "basis" in
+let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
+
+let charge =
+  match Command_line.get "charge" with
+  | Some x -> ( if x.[0] = 'm' then
+                  ~- (int_of_string (String.sub x 1 (String.length x - 1)))
+                else
+                  int_of_string x )
+  | None   -> 0
+in
+    
+let multiplicity =
+  match Command_line.get "multiplicity" with
+  | Some x -> int_of_string x
+  | None -> 1
+in
+   #+END_SRC
+
+* Computation
+  We first read the =xyz= file to create a molecule:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.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:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml  
+let ao_basis =
+  Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file
+in
+  #+END_SRC
+
+  We create a simulation from the nuclei and the basis set:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml  
+let simulation = Qcaml.Simulation.make ~nuclei ao_basis in
+  #+END_SRC
+
+  and we can make the Hartree-Fock computation:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml  
+let hf = Qcaml.Mo.Hartree_fock.make ~guess:`Huckel simulation in
+  #+END_SRC
+
+* Output
+
+  We print the convergence of the calculation:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_hartree_fock.ml  
+Format.printf "@[%a@]" (Mo.Hartree_fock.pp) hf
+  #+END_SRC
+
+  
diff --git a/examples/ex_integrals.org b/examples/ex_integrals.org
index 0ce91cf..241c3b6 100644
--- a/examples/ex_integrals.org
+++ b/examples/ex_integrals.org
@@ -3,7 +3,7 @@
 #+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
+molecule 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.
 
diff --git a/linear_algebra/lib/matrix.ml b/linear_algebra/lib/matrix.ml
index f01f700..b30c853 100644
--- a/linear_algebra/lib/matrix.ml
+++ b/linear_algebra/lib/matrix.ml
@@ -176,17 +176,17 @@ let gemm ?m ?n ?k ?beta ?c ?(transa=`N) ?(alpha=1.0) a ?(transb=`N) b =
   in
   gemm ?m ?n ?k ?beta ?c ~transa ~alpha a ~transb b
 
-let gemm_nn ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b =
-  gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`N ~transb:`N
+let gemm_nn ?m ?n ?k ?beta ?c ?(alpha=1.0) a b =
+  gemm ?m ?n ?k ?beta ?c ~alpha a b ~transa:`N ~transb:`N
 
-let gemm_nt ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b =
-  gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`N ~transb:`T
+let gemm_nt ?m ?n ?k ?beta ?c ?(alpha=1.0) a b =
+  gemm ?m ?n ?k ?beta ?c ~alpha a b ~transa:`N ~transb:`T
 
-let gemm_tn ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b =
-  gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`T ~transb:`N
+let gemm_tn ?m ?n ?k ?beta ?c ?(alpha=1.0) a b =
+  gemm ?m ?n ?k ?beta ?c ~alpha a b ~transa:`T ~transb:`N
 
-let gemm_tt ?m ?n ?k ?(beta=0.) ?c ?(alpha=1.0) a b =
-  gemm ?m ?n ?k ~beta ?c ~alpha a b ~transa:`T ~transb:`T
+let gemm_tt ?m ?n ?k ?beta ?c ?(alpha=1.0) a b =
+  gemm ?m ?n ?k ?beta ?c ~alpha a b ~transa:`T ~transb:`T
 
 
 let gemm_trace ?(transa=`N) a ?(transb=`N) b =
diff --git a/mo/test/hf.ml b/mo/test/hf.ml
new file mode 100644
index 0000000..72925b0
--- /dev/null
+++ b/mo/test/hf.ml
@@ -0,0 +1,24 @@
+open Alcotest
+open Particles
+    
+let wd = Common.Qcaml.root ^ Filename.dir_sep ^ "test"
+
+let tests = 
+  [ "HF Water", `Quick, fun () -> 
+    let nuclei =
+      wd ^ Filename.dir_sep ^ "water.xyz"
+      |> Nuclei.of_xyz_file 
+    in
+    let basis_filename =
+      wd ^ Filename.dir_sep ^ "cc-pvdz"
+    in
+    let ao_basis = 
+      Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian
+        ~cartesian:false ~nuclei basis_filename
+    in
+    let simulation = Simulation.make ~nuclei ao_basis in
+    let hf = Mo.Hartree_fock.make ~guess:`Huckel simulation in
+    Format.printf "%a" (Mo.Hartree_fock.pp) hf;
+    check (float 1.e-10) "Energy" (-76.0267987006) (Mo.Hartree_fock.energy hf);
+  ]
+    
diff --git a/qcaml/lib/dune b/qcaml/lib/dune
index 9453e6b..0656da0 100644
--- a/qcaml/lib/dune
+++ b/qcaml/lib/dune
@@ -4,11 +4,13 @@
  (name qcaml)
  (public_name qcaml)
  (libraries
+	qcaml.ao
         qcaml.common
-        qcaml.particles
         qcaml.gaussian
         qcaml.gaussian_integrals
+	qcaml.mo
 	qcaml.operators
-	qcaml.ao
+        qcaml.particles
+	qcaml.simulation
  )
  (synopsis "Main QCaml entry point"))
diff --git a/qcaml/lib/qcaml.ml b/qcaml/lib/qcaml.ml
index 8784bd8..13dcd04 100644
--- a/qcaml/lib/qcaml.ml
+++ b/qcaml/lib/qcaml.ml
@@ -1,8 +1,9 @@
-module Common = Common
 module Ao = Ao
+module Common = Common
 module Gaussian = Gaussian
 module Gaussian_integrals = Gaussian_integrals
 module Linear_algebra = Linear_algebra
+module Mo = Mo
 module Operators = Operators
 module Particles = Particles
-
+module Simulation = Simulation
diff --git a/test/run_tests.ml b/test/run_tests.ml
index fecebd5..55b8f39 100644
--- a/test/run_tests.ml
+++ b/test/run_tests.ml
@@ -12,6 +12,7 @@ let test_suites: unit Alcotest.test list = [
   "Ao_basis.Ao_basis_gaussian", Test_ao_basis.Ao_basis_gaussian.tests;
   "Ao_basis.Ao_basis", Test_ao_basis.Ao_basis.tests;
   "Mo_basis.Guess", Test_mo_basis.Guess.tests;
+  "Mo_basis.Hartree_Fock", Test_mo_basis.Hf.tests;
 ]
 
 let () = Alcotest.run "QCaml" test_suites