diff --git a/Basis/BasisLexer.mll b/Basis/BasisLexer.mll
new file mode 100644
index 0000000..9162984
--- /dev/null
+++ b/Basis/BasisLexer.mll
@@ -0,0 +1,33 @@
+{
+exception SyntaxError of string
+
+open GamessParser
+
+}
+
+let eol = ['\n']
+let white = [' ' '\t']+
+let element   = ['A'-'Z' 'a'-'z']+ white? eol
+let ang_mom  = ['S' 'P' 'D' 'F' 'G' 'H' 'I' 'J' 'K' 'L' 'M' 'N' 'O' 
+                's' 'p' 'd' 'f' 'g' 'h' 'i' 'j' 'k' 'l' 'm' 'n' 'o' ]
+               white
+let integer  = ['0'-'9']+ 
+let real = '-'? integer '.' integer (['e' 'E'] ('+'|'-')? integer)? 
+
+
+rule read_all_rule  = parse
+  | eol            { EOL }
+  | white          { SPACE }
+  | ang_mom as a   { ANG_MOM (a.[0]) }
+  | element as e   { ELEMENT (String.trim e) }
+  | integer as i   { INTEGER (int_of_string i) }
+  | real    as f   { FLOAT (float_of_string f) }
+  | eof            { EOF }
+
+
+{
+  let rec read_all lexbuf =
+     match read_all_rule lexbuf with
+     | SPACE -> read_all_rule lexbuf
+     | x -> x
+}
diff --git a/Basis/ContractedShell.ml b/Basis/ContractedShell.ml
new file mode 100644
index 0000000..47ebe42
--- /dev/null
+++ b/Basis/ContractedShell.ml
@@ -0,0 +1,127 @@
+open Util
+open Constants
+open Coordinate
+
+type t = {
+  expo      : float array;         (* Gaussian exponents *)
+  coef      : float array;         (* Contraction coefficients *)
+  center    : Coordinate.t;        (* Center of all the Gaussians *)
+  totAngMom : AngularMomentum.t;  (* Total angular momentum *)
+  size      : int;                 (* Number of contracted Gaussians *)
+  norm_coef : float array;         (* Normalization coefficient of the class
+                                      corresponding to the i-th contraction *)
+  norm_coef_scale : float array;   (* Inside a class, the norm is the norm
+                                      of the function with (totAngMom,0,0) *. 
+                                      this scaling factor *)
+  index    : int;                  (* Index in the array of contracted shells *)
+  powers   : Zkey.t array;         (* Array of Zkeys corresponding to the
+                                      powers of (x,y,z) in the class *)
+}
+
+module Am = AngularMomentum
+
+(** Normalization coefficient of contracted function i, which depends on the
+    exponent and the angular momentum. Two conventions can be chosen : a single
+    normalisation factor for all functions of the class, or a coefficient which
+    depends on the powers of x,y and z.
+    Returns, for each contracted function, an array of functions taking as
+    argument the [|x;y;z|] powers.
+*)
+let compute_norm_coef expo totAngMom =
+  let atot =
+    Am.to_int totAngMom
+  in
+  let factor int_array = 
+       let dfa = Array.map (fun j ->
+           ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) 
+         ) int_array
+       in
+       sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
+  in
+  let expo = 
+    if atot mod 2 = 0 then
+      Array.map (fun alpha ->
+        let alpha_2 = alpha +. alpha in
+        (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
+      ) expo
+    else
+      Array.map (fun alpha ->
+        let alpha_2 = alpha +. alpha in
+        (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
+      ) expo
+  in
+  Array.map (fun x -> let f a = x *. (factor a) in f) expo
+
+
+let make ~index ~expo ~coef ~center ~totAngMom =
+  assert (Array.length expo = Array.length coef);
+  assert (Array.length expo > 0);
+  let norm_coef_func =
+    compute_norm_coef expo totAngMom
+  in
+  let powers =
+    Am.zkey_array (Am.Singlet totAngMom) 
+  in
+  let norm_coef =
+    Array.map (fun f -> f [| Am.to_int totAngMom ; 0 ; 0 |])  norm_coef_func
+  in
+  let norm_coef_scale = 
+    Array.map (fun a ->
+      (norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0)
+    ) powers
+  in
+  { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
+    norm_coef_scale ; powers }
+
+
+let with_index a i =
+  { a with index = i }
+
+
+let to_string s =
+  let coord = s.center in
+  let open Printf in
+  (match s.totAngMom with
+   | Am.S -> sprintf "%3d    " (s.index+1)
+   | _                  -> sprintf "%3d-%-3d" (s.index+1) (s.index+(Array.length s.powers))
+  ) ^
+  ( sprintf "%1s %8.3f %8.3f %8.3f " (Am.to_string s.totAngMom)
+      (get X coord) (get Y coord) (get Z coord) ) ^
+  (Array.map2 (fun e c -> sprintf "%16.8e  %16.8e" e c) s.expo s.coef
+   |> Array.to_list |> String.concat (sprintf "\n%36s" " ") ) 
+
+
+(** Normalization coefficient of contracted function i, which depends on the
+    exponent and the angular momentum. Two conventions can be chosen : a single
+    normalisation factor for all functions of the class, or a coefficient which
+    depends on the powers of x,y and z.
+    Returns, for each contracted function, an array of functions taking as
+    argument the [|x;y;z|] powers.
+*)
+let compute_norm_coef expo totAngMom =
+  let atot =
+    Am.to_int totAngMom
+  in
+  let factor int_array = 
+       let dfa = Array.map (fun j ->
+           (float_of_int (1 lsl j) *. fact j) /. fact (j+j) 
+         ) int_array
+       in
+       sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
+  in
+  let expo = 
+    if atot mod 2 = 0 then
+      Array.map (fun alpha ->
+        let alpha_2 = alpha +. alpha in
+        (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
+      ) expo
+    else
+      Array.map (fun alpha ->
+        let alpha_2 = alpha +. alpha in
+        (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
+      ) expo
+  in
+  Array.map (fun x -> let f a = x *. factor a in f) expo
+
+
+
diff --git a/Basis/ContractedShell.mli b/Basis/ContractedShell.mli
new file mode 100644
index 0000000..714fae7
--- /dev/null
+++ b/Basis/ContractedShell.mli
@@ -0,0 +1,26 @@
+type t = private {
+  expo      : float array;
+  coef      : float array;
+  center    : Coordinate.t;
+  totAngMom : AngularMomentum.t;
+  size      : int;
+  norm_coef : float array;
+  norm_coef_scale : float array;
+  index    : int;
+  powers   : Zkey.t array;
+}
+
+
+
+(** Pretty-printing of the contracted shell in a string *)
+val to_string : t -> string
+
+(** Creates a contracted shell *)
+val make :
+  index:int ->
+  expo:float array ->
+  coef:float array ->
+  center:Coordinate.t -> totAngMom:AngularMomentum.t -> t
+
+(** Returns a copy of the contracted shell with a modified index *)
+val with_index  : t -> int -> t
diff --git a/Basis/GamessParser.mly b/Basis/GamessParser.mly
new file mode 100644
index 0000000..35e3c68
--- /dev/null
+++ b/Basis/GamessParser.mly
@@ -0,0 +1,49 @@
+/* Parses basis sets GAMESS format */
+
+%{
+
+%}
+
+%token <string> ELEMENT
+%token <char> ANG_MOM
+%token <int> INTEGER
+%token <float> FLOAT
+%token SPACE
+%token EOL
+%token EOF
+
+%start input
+%type <GeneralBasis.t> input
+
+%% /* Grammar rules and actions follow */
+
+input:
+  | basis { $1 }
+  | EOL input { $2 }
+
+basis:
+  | element shell_array EOL { ($1, $2) }
+  | element shell_array EOF { ($1, $2) }
+
+element:
+  | ELEMENT { Element.of_string $1 }
+
+shell_array:
+  | shell_list { Array.of_list @@ List.rev $1 }
+
+shell_list:
+  | { [] }
+  | shell_list shell { $2 :: $1 }
+
+shell:
+  | ANG_MOM INTEGER EOL primitive_list { (AngularMomentum.of_char $1, Array.of_list @@ List.rev $4 ) } 
+
+primitive_list: 
+  | { [] }
+  | primitive_list primitive { $2 :: $1 }
+
+primitive: 
+  | INTEGER FLOAT FLOAT EOL { GeneralBasis.{exponent=$2 ; coefficient=$3 } }
+
+
+
diff --git a/Basis/GamessReader.ml b/Basis/GamessReader.ml
new file mode 100644
index 0000000..eca4fb9
--- /dev/null
+++ b/Basis/GamessReader.ml
@@ -0,0 +1,18 @@
+(** Read a basis set file in GAMESS format and return an association list where the key is an
+    Element.t and the value is the parsed basis set. *)
+let read ~filename =
+  let lexbuf =
+    let ic = open_in filename in
+    Lexing.from_channel ic
+  in
+  let rec aux accu = 
+    try
+      let key, basis =
+        GamessParser.input BasisLexer.read_all lexbuf
+      in
+      aux ((key, basis)::accu)
+    with
+    | Parsing.Parse_error -> List.rev accu
+  in
+  aux []
+
diff --git a/Basis/GeneralBasis.ml b/Basis/GeneralBasis.ml
new file mode 100644
index 0000000..470c3d6
--- /dev/null
+++ b/Basis/GeneralBasis.ml
@@ -0,0 +1,40 @@
+(** General basis set read from a file  *)
+type primitive = {
+    exponent: float ;
+    coefficient: float
+    }
+
+type general_contracted_shell = AngularMomentum.t * (primitive array)
+
+type t = Element.t * (general_contracted_shell array)
+
+
+module Am = AngularMomentum
+
+let string_of_primitive ?id prim =
+  match id with
+  | None -> (string_of_float prim.exponent)^" "^(string_of_float prim.coefficient)
+  | Some i -> (string_of_int i)^" "^(string_of_float prim.exponent)^" "^(string_of_float prim.coefficient)
+
+
+let string_of_contracted_shell (angular_momentum, prim_array) =
+  let n = 
+    Array.length prim_array
+  in
+  Printf.sprintf "%s %d\n%s"
+    (Am.to_string angular_momentum) n
+    (Array.init n (fun i -> string_of_primitive ~id:(i+1) prim_array.(i))
+     |> Array.to_list
+     |> String.concat "\n")
+
+
+let string_of_contracted_shell_array a =
+  Array.map string_of_contracted_shell a
+  |> Array.to_list
+  |> String.concat "\n"
+
+
+let to_string (name, contracted_shell_array) =
+  Printf.sprintf "%s\n%s" name (string_of_contracted_shell_array contracted_shell_array)
+
+
diff --git a/Basis/Orthonormalization.ml b/Basis/Orthonormalization.ml
new file mode 100644
index 0000000..82b2607
--- /dev/null
+++ b/Basis/Orthonormalization.ml
@@ -0,0 +1,29 @@
+open Util
+open Lacaml.D
+
+type t = 
+| Lowdin    of Mat.t
+| Canonical of Mat.t
+| Svd       of Mat.t
+
+
+let make_lowdin ?(thresh=1.e-12) ~overlap =
+  
+  let u_vec, u_val = diagonalize_symm overlap in
+
+  Vec.iter (fun x -> if x < thresh then
+      invalid_arg "Orthonormalization.make_lowdin") u_val;
+
+  let u_val = Vec.reci (Vec.sqrt u_val) in
+
+  let u_vec' =
+    Mat.init_cols (Mat.dim1 u_vec) (Mat.dim2 u_vec) (fun i j -> u_vec.{i,j} *. u_val.{j})
+  in
+  let result = 
+    gemm u_vec' ~transb:`T u_vec
+  in
+
+  Lowdin result
+     
+
+let make = make_lowdin
diff --git a/HartreeFock/Fock.ml b/HartreeFock/Fock.ml
new file mode 100644
index 0000000..94b4481
--- /dev/null
+++ b/HartreeFock/Fock.ml
@@ -0,0 +1,33 @@
+open Lacaml.D
+open Simulation
+
+type t = Mat.t
+
+let make ~density simulation =
+  let m_P  = density
+  and m_T  = Lazy.force simulation.kin_ints
+  and m_V  = Lazy.force simulation.eN_ints
+  and m_G  = Lazy.force simulation.ee_ints
+  in
+  let nBas = Mat.dim1 m_T
+  in
+
+  let m_F = Mat.add m_T m_V in
+  for sigma = 1 to nBas do
+    for nu = 1 to nBas do
+      for lambda = 1 to nBas do
+        let p = m_P.{lambda,sigma} in
+        for mu = 1 to nu do
+          m_F.{mu,nu} <- m_F.{mu,nu} +.  p *. 
+                        (m_G.{mu,lambda,nu,sigma} -. 0.5 *. m_G.{mu,lambda,sigma,nu})
+        done
+      done
+    done
+  done;
+  for nu = 1 to nBas do
+    for mu = 1 to nu do
+      m_F.{nu,mu} <- m_F.{mu,nu}
+    done
+  done;
+  m_F
+
diff --git a/HartreeFock/Guess.ml b/HartreeFock/Guess.ml
new file mode 100644
index 0000000..8a4c471
--- /dev/null
+++ b/HartreeFock/Guess.ml
@@ -0,0 +1,17 @@
+open Lacaml.D
+
+type guess = 
+| Hcore of Mat.t
+
+type t = guess
+
+
+let make ?guess:(guess=`Hcore) simulation = 
+  let eN_ints  = Lazy.force simulation.Simulation.eN_ints
+  and kin_ints = Lazy.force simulation.Simulation.kin_ints
+  in
+  match guess with
+  | `Hcore -> Hcore (Mat.add eN_ints kin_ints)
+
+
+
diff --git a/HartreeFock/Guess.mli b/HartreeFock/Guess.mli
new file mode 100644
index 0000000..499822d
--- /dev/null
+++ b/HartreeFock/Guess.mli
@@ -0,0 +1,12 @@
+(** Guess for Hartree-Fock calculations. *)
+
+type guess = 
+| Hcore of Lacaml.D.Mat.t
+
+type t = guess
+
+
+val make : ?guess:[ `Hcore ] -> Simulation.t -> t
+
+
+
diff --git a/HartreeFock/HartreeFock.ml b/HartreeFock/HartreeFock.ml
new file mode 100644
index 0000000..4673eac
--- /dev/null
+++ b/HartreeFock/HartreeFock.ml
@@ -0,0 +1,9 @@
+open Util
+open Simulation
+
+let make ?guess simulation =
+  if simulation.electrons.Electrons.multiplicity = 1 then
+    RHF.make ?guess simulation
+  else
+    invalid_arg "UHF or ROHF not implemented"
+
diff --git a/HartreeFock/HartreeFock_type.ml b/HartreeFock/HartreeFock_type.ml
new file mode 100644
index 0000000..0d60f72
--- /dev/null
+++ b/HartreeFock/HartreeFock_type.ml
@@ -0,0 +1,11 @@
+type t =
+  {
+    guess        : Guess.t;
+    eigenvectors : Lacaml.D.Mat.t ;
+    eigenvalues  : Lacaml.D.Vec.t ;
+    energy       : float ;
+    iterations   : float array;
+  }
+
+
+
diff --git a/HartreeFock/HartreeFock_type.mli b/HartreeFock/HartreeFock_type.mli
new file mode 100644
index 0000000..7cdbcaa
--- /dev/null
+++ b/HartreeFock/HartreeFock_type.mli
@@ -0,0 +1,10 @@
+type t = 
+  {
+    guess        : Guess.t;
+    eigenvectors : Lacaml.D.Mat.t ;
+    eigenvalues  : Lacaml.D.Vec.t ;
+    energy       : float ;
+    iterations   : float array;
+  }
+
+
diff --git a/HartreeFock/RHF.ml b/HartreeFock/RHF.ml
new file mode 100644
index 0000000..64df9e3
--- /dev/null
+++ b/HartreeFock/RHF.ml
@@ -0,0 +1,126 @@
+open Util
+open Lacaml.D
+open Simulation
+
+let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
+  ?threshold_SCF:(threshold_SCF=1.e-6) simulation =
+
+  (* Number of occupied MOs *)
+  let nocc =
+    simulation.electrons.Electrons.n_alpha
+  in
+
+  (* Initial guess *)
+  let guess = 
+    Guess.make ~guess simulation 
+  in
+
+  (* Orthogonalization matrix *)
+  let m_X = 
+    match Lazy.force simulation.overlap_ortho with
+    | Lowdin    x -> x
+    | Canonical x -> x
+    | Svd       x -> x
+  in
+
+  (* Core Hamiltonian *)
+  let m_Hc =
+    let m_T = Lazy.force simulation.kin_ints
+    and m_V = Lazy.force simulation.eN_ints
+    in Mat.add m_T m_V
+  in
+
+  (* Overlap matrix *)
+  let m_S =
+    Lazy.force  simulation.overlap
+  in
+
+
+  (* SCF iterations *)
+  let rec loop nSCF iterations m_C =
+
+      (* Density matrix over nocc occupied MOs *)
+      let m_P =
+        gemm  ~alpha:2.  ~transb:`T  ~k:nocc  m_C  m_C
+      in
+
+      (* Fock matrix in AO basis *)
+      let m_F =
+        Fock.make  ~density:m_P  simulation
+      in
+
+      (* Fock matrix in MO basis *)
+      let m_Fmo =
+        xt_o_x  m_F  m_X
+      in
+
+      (* MOs in old MO basis *)
+      let m_C', eigenvalues =
+        diagonalize_symm  m_Fmo
+      in
+
+      (* MOs in AO basis *)
+      let m_C =
+        gemm  m_X  m_C'
+      in
+
+      (* Hartree-Fock energy *)
+      let energy =
+        simulation.nuclear_repulsion +. 0.5 *.
+          Mat.gemm_trace  m_P  (Mat.add  m_Hc  m_F)
+      in
+
+      (* Convergence criterion *)
+      let commutator =
+        let fps =
+          gemm  m_F  (gemm  m_P  m_S)
+        and spf =
+          gemm  m_S  (gemm  m_P  m_F)
+        in
+        Mat.sub  fps  spf
+        |> Mat.as_vec
+        |> amax
+      in
+
+      let converged =
+        nSCF = max_scf  ||  (abs_float commutator) < threshold_SCF
+      in
+
+      Printf.printf "%d  %16.10f  %10.4e\n%!" nSCF energy commutator;
+
+      if not converged then
+          loop (nSCF+1) (energy :: iterations) m_C
+      else
+          let iterations = 
+            List.rev iterations
+            |> Array.of_list
+          in
+          { HartreeFock_type.
+            guess ;
+            eigenvectors = m_C ;
+            eigenvalues ;
+            energy = iterations.(Array.length iterations - 1) ;
+            iterations ;
+          }
+  in
+
+
+  (* Guess coefficients *)
+  let m_H = 
+    match guess with
+    | Guess.Hcore m_H -> m_H 
+  in
+  let m_Hmo =
+    xt_o_x  m_H  m_X
+  in
+  let m_C', _ =
+    diagonalize_symm  m_Hmo
+  in
+  let m_C =
+    gemm  m_X  m_C'
+  in
+
+  loop 1 [] m_C
+
+
+
diff --git a/Utils/AngularMomentum.ml b/Utils/AngularMomentum.ml
new file mode 100644
index 0000000..d7877a2
--- /dev/null
+++ b/Utils/AngularMomentum.ml
@@ -0,0 +1,134 @@
+open Powers
+
+exception AngularMomentumError of string
+
+type t =
+  | S | P | D | F | G | H | I | J | K | L | M | N | O
+
+let of_char = function
+  | 's' | 'S' -> S      | 'p' | 'P' -> P
+  | 'd' | 'D' -> D      | 'f' | 'F' -> F
+  | 'g' | 'G' -> G      | 'h' | 'H' -> H
+  | 'i' | 'I' -> I      | 'j' | 'J' -> J
+  | 'k' | 'K' -> K      | 'l' | 'L' -> L
+  | 'm' | 'M' -> M      | 'n' | 'N' -> N
+  | 'o' | 'O' -> O
+  | c -> raise (AngularMomentumError (String.make 1 c))
+
+let to_string = function
+  | S -> "S"       | P -> "P"
+  | D -> "D"       | F -> "F"
+  | G -> "G"       | H -> "H"
+  | I -> "I"       | J -> "J"
+  | K -> "K"       | L -> "L"
+  | M -> "M"       | N -> "N"
+  | O -> "O"
+
+let to_char = function
+  | S -> 'S'       | P -> 'P'
+  | D -> 'D'       | F -> 'F'
+  | G -> 'G'       | H -> 'H'
+  | I -> 'I'       | J -> 'J'
+  | K -> 'K'       | L -> 'L'
+  | M -> 'M'       | N -> 'N'
+  | O -> 'O'
+
+let to_int = function
+  | S -> 0         | P -> 1
+  | D -> 2         | F -> 3
+  | G -> 4         | H -> 5
+  | I -> 6         | J -> 7
+  | K -> 8         | L -> 9
+  | M -> 10        | N -> 11
+  | O -> 12
+
+let of_int = function
+  | 0 -> S         | 1 -> P
+  | 2 -> D         | 3 -> F
+  | 4 -> G         | 5 -> H
+  | 6 -> I         | 7 -> J
+  | 8 -> K         | 9 -> L
+  | 10 -> M        | 11 -> N
+  | 12 -> O
+  | c -> raise (AngularMomentumError (string_of_int c))
+
+
+
+type kind =
+| Singlet of  t
+| Doublet of (t*t)
+| Triplet of (t*t*t)
+| Quartet of (t*t*t*t)
+
+
+let n_functions a =
+  let a = 
+    to_int a
+  in
+  (a*a + 3*a + 2)/2
+
+
+(** Returns an array of Zkeys corresponding to all possible angular momenta *)
+let zkey_array a = 
+  let keys_1d l =
+    let create_z { x ; y ; _ } =
+      Powers.of_int_tuple (x,y,l-(x+y))
+    in
+    let rec create_y accu xyz =
+      let { x ; y ; z } = xyz in
+      match y with
+      | 0 -> (create_z xyz)::accu
+      | i -> let ynew = y-1 in
+        create_y ( (create_z xyz)::accu) (Powers.of_int_tuple (x,ynew,z))
+    in
+    let rec create_x accu xyz =
+      let { x ; y ; z } = xyz in
+      match x with
+      | 0 -> (create_y [] xyz)@accu
+      | i -> let xnew = x-1 in
+        let ynew = l-xnew in
+        create_x ((create_y [] xyz)@accu) (Powers.of_int_tuple (xnew, ynew, z))
+    in
+    create_x [] (Powers.of_int_tuple (l,0,0))
+    |> List.rev
+  in
+
+  begin
+    match a with
+    | Singlet l1 -> 
+      List.map (fun x -> Zkey.of_powers (Zkey.Three x)) (keys_1d @@ to_int l1)
+
+    | Doublet (l1, l2) -> 
+      List.map (fun a ->
+          List.map (fun b ->
+              Zkey.of_powers (Zkey.Six (a,b))) (keys_1d @@ to_int l2)
+        ) (keys_1d @@ to_int l1)
+      |> List.concat
+
+    | Triplet (l1, l2, l3) ->
+
+      List.map (fun a ->
+          List.map (fun b ->
+              List.map (fun c ->
+                  Zkey.of_powers (Zkey.Nine (a,b,c))) (keys_1d @@ to_int l3)
+            ) (keys_1d @@ to_int l2)
+          |> List.concat
+        ) (keys_1d @@ to_int l1)
+      |> List.concat
+
+    | Quartet (l1, l2, l3, l4) ->
+
+      List.map (fun a ->
+          List.map (fun b ->
+              List.map (fun c ->
+                  List.map (fun d ->
+                      Zkey.of_powers (Zkey.Twelve (a,b,c,d))) (keys_1d @@ to_int l4)
+                ) (keys_1d @@ to_int l3)
+              |> List.concat
+            ) (keys_1d @@ to_int l2)
+          |> List.concat
+        ) (keys_1d @@ to_int l1)
+        |> List.concat
+  end
+  |> Array.of_list
+
diff --git a/Utils/Electrons.ml b/Utils/Electrons.ml
new file mode 100644
index 0000000..b951d91
--- /dev/null
+++ b/Utils/Electrons.ml
@@ -0,0 +1,23 @@
+(** Number of alpha and beta electrons *)
+
+type t = {
+  n_alpha : int ;
+  n_beta  : int ;
+  multiplicity : int;
+}
+
+
+let make ?multiplicity:(multiplicity=1) ?charge:(charge=0) nuclei =
+  let positive_charges =
+    Array.fold_left (fun accu (e, _) -> accu + Charge.to_int (Element.to_charge e) )
+      0 nuclei
+  in
+  let negative_charges = charge - positive_charges in
+  let n_elec = - negative_charges in
+  let n_beta = ((n_elec - multiplicity)+1)/2 in
+  let n_alpha = n_elec - n_beta in
+  let result = { n_alpha ; n_beta ; multiplicity } in
+  if multiplicity <> (n_alpha - n_beta)+1 then
+    invalid_arg "Electrons.make";
+  result
+
diff --git a/Utils/PositiveFloat.ml b/Utils/PositiveFloat.ml
new file mode 100644
index 0000000..fbf46c0
--- /dev/null
+++ b/Utils/PositiveFloat.ml
@@ -0,0 +1,14 @@
+type t = float
+
+let of_float x =
+  assert ( x >= 0. );
+  x
+
+external to_float : t -> float = "%identity"
+
+let to_string x =
+  let f = to_float x in string_of_float f
+
+let of_string x =
+ let f = float_of_string x in of_float f
+
diff --git a/Utils/PositiveFloat.mli b/Utils/PositiveFloat.mli
new file mode 100644
index 0000000..2d3d60b
--- /dev/null
+++ b/Utils/PositiveFloat.mli
@@ -0,0 +1,5 @@
+type t = private float
+val of_float : float -> t
+val to_float : t -> float
+val to_string : t -> string
+val of_string : string -> t
diff --git a/run_hartree_fock.ml b/run_hartree_fock.ml
new file mode 100644
index 0000000..144da8d
--- /dev/null
+++ b/run_hartree_fock.ml
@@ -0,0 +1,53 @@
+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 charge      : int ref = ref 0
+let multiplicity: int ref = ref 1
+
+
+let speclist = [
+  ( "-b" , Arg.String (fun x -> basis_file := Some x),
+    "File containing the atomic basis set") ;
+  ( "-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") ;
+  ( "-charge" , Arg.Int (fun x -> charge := x),
+    "Charge of the system") ;
+  ( "-mult" , Arg.Int (fun x -> multiplicity := x),
+    "Spin multiplicity of the system") ;
+]
+
+let run ~out =
+  (*
+  let gc = Gc.get () in
+  Gc.set { gc with minor_heap_size=(262144 / 16) };
+  *)
+  let 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
+  and charge = !charge
+  and multiplicity = !multiplicity
+  in
+
+  let s =
+    Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file 
+  in
+
+  let _ =
+    HartreeFock.make s 
+  in
+  Printf.printf "Done.\n%!";
+  ()
+
+
+let () =
+  let usage_msg = "Available options:" in
+  Arg.parse speclist (fun _ -> ()) usage_msg;
+  run ~out:!out_file
+