From ff872e94794126e507c18354d1b181c75d9171f6 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 28 Jun 2023 16:22:56 +0200
Subject: [PATCH] Slower

---
 common/charge.org                             |  96 -----
 common/command_line.org                       | 354 ------------------
 common/constants.org                          |  79 ----
 common/lib/charge.ml                          |  13 -
 common/lib/charge.mli                         |  22 +-
 common/lib/command_line.ml                    | 100 +----
 common/lib/command_line.mli                   |  94 ++++-
 common/lib/constants.ml                       |  28 +-
 common/lib/constants.mli                      |  41 +-
 common/lib/smallarray.ml                      | 103 +++++
 common/lib/smallarray.mli                     |  24 ++
 .../lib/two_electron_rr_vectorized.ml         | 217 ++++-------
 top/lib/install_printers.ml                   |   1 -
 13 files changed, 344 insertions(+), 828 deletions(-)
 delete mode 100644 common/charge.org
 delete mode 100644 common/command_line.org
 delete mode 100644 common/constants.org
 create mode 100644 common/lib/smallarray.ml
 create mode 100644 common/lib/smallarray.mli

diff --git a/common/charge.org b/common/charge.org
deleted file mode 100644
index cef7097..0000000
--- a/common/charge.org
+++ /dev/null
@@ -1,96 +0,0 @@
-#+begin_src elisp tangle: no :results none :exports none
-(setq pwd (file-name-directory buffer-file-name))
-(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
-(setq lib (concat pwd "lib/"))
-(setq testdir (concat pwd "test/"))
-(setq mli (concat lib name ".mli"))
-(setq ml  (concat lib name ".ml"))
-(setq test-ml  (concat testdir name ".ml"))
-(org-babel-tangle)
-#+end_src 
-
-
-* Charge
-  :PROPERTIES:
-  :header-args: :noweb yes :comments both
-  :END:
-
-** Type
-
-   <<<~Charge.t~>>>
-   #+begin_src ocaml :tangle (eval mli)
-type t
-   #+end_src
-
-   This type should be used for all charges in the program (electrons, nuclei,...).
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-type t = float
-   #+end_src
-
-** Conversions
-
-   #+begin_src ocaml :tangle (eval mli)
-val of_float : float -> t
-val to_float : t -> float
-
-val of_int   : int -> t
-val to_int   : t -> int
-
-val of_string: string -> t
-val to_string: t -> string 
-   #+end_src
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-external of_float : float -> t = "%identity"
-external to_float : t -> float = "%identity"
-
-let of_int = float_of_int
-let to_int = int_of_float
-
-let of_string = float_of_string 
-
-let to_string x = 
-  if x >  0. then
-    Printf.sprintf "+%f" x
-  else if x < 0. then
-    Printf.sprintf "%f" x
-  else 
-    "0.0"
-   #+end_src
-
-** Simple operations
-
-   #+begin_src ocaml :tangle (eval mli)
-val ( + ) : t -> t -> t
-val ( - ) : t -> t -> t
-val ( * ) : t -> float -> t
-val ( / ) : t -> float -> t 
-val is_null : t -> bool
-   #+end_src
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let gen_op op =
-  fun a b ->
-  op (to_float a) (to_float b)
-  |> of_float
-
-let ( + ) = gen_op ( +. )
-let ( - ) = gen_op ( -. )
-let ( * ) = gen_op ( *. )
-let ( / ) = gen_op ( /. )
-
-let is_null t = t == 0.
-   #+end_src
-
-** Printers
-
-   #+begin_src ocaml :tangle (eval mli)
-val pp : Format.formatter -> t -> unit
-   #+end_src
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let pp ppf x = 
-  Format.fprintf ppf "@[%s@]" (to_string x)
-   #+end_src
-
diff --git a/common/command_line.org b/common/command_line.org
deleted file mode 100644
index 17ac72e..0000000
--- a/common/command_line.org
+++ /dev/null
@@ -1,354 +0,0 @@
-#+begin_src elisp tangle: no :results none :exports none
-(setq pwd (file-name-directory buffer-file-name))
-(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
-(setq lib (concat pwd "lib/"))
-(setq testdir (concat pwd "test/"))
-(setq mli (concat lib name ".mli"))
-(setq ml  (concat lib name ".ml"))
-(setq test-ml  (concat testdir name ".ml"))
-(org-babel-tangle)
-#+end_src
-
-* Command line
-  :PROPERTIES:
-  :header-args: :noweb yes :comments both
-  :END:
-
-  This module is a wrapper around the ~Getopt~ library and helps to
-  define command-line arguments.
-
-  Here is an example of how to use this module.
-  First, define the specification:
-  #+begin_src ocaml :tangle no
-let open Command_line in
-begin
-  set_header_doc (Sys.argv.(0) ^ " - One-line description");
-  set_description_doc "Long description of the command.";
-  set_specs
-    [ { short='c'; long="check"; opt=Optional;
-        doc="Checks the input data";
-        arg=Without_arg; };
-
-      { short='b' ; long="basis" ; opt=Mandatory;
-        arg=With_arg "<string>";
-        doc="Name of the file containing the basis set"; } ;
-
-      { short='m' ; long="multiplicity" ; opt=Optional;
-        arg=With_arg "<int>";
-        doc="Spin multiplicity (2S+1). Default is singlet"; } ;
-    ]
-end;
-  #+end_src
-
-  Then, define what to do with the arguments:
-  #+begin_src ocaml :tangle no
-let c =
-  Command_line.get_bool "check"
-in
-
-let basis =
-  match Command_line.get "basis" with
-  | Some x -> x
-  | None -> assert false
-in
-
-let multiplicity =
-  match Command_line.get "multiplicity" with
-  | None -> 1
-  | Some n -> int_of_string n
-in
-  #+end_src
-
-** Types
-
-   #+NAME:type
-   #+begin_src ocaml :tangle (eval mli)
-type short_opt     = char
-type long_opt      = string
-type optional      = Mandatory | Optional
-type documentation = string
-type argument      = | With_arg of string
-                     | Without_arg
-                     | With_opt_arg of string
-
-type description = {
-  short: short_opt ;
-  long : long_opt  ;
-  opt  : optional  ;
-  doc  : documentation ;
-  arg  : argument ;
-}
-
-   #+end_src
-
-   - <<<Short option>>>: in the command line, a dash with a single character
-     (ex: =ls -l=)
-   - <<<Long option>>>: in the command line, two dashes with a word
-     (ex: =ls --directory=)
-   - Command-line options can be ~Mandatory~ or ~Optional~
-   - Documentation of the option is used in the help function
-   - Some options require an argument (~ls --ignore="*.ml"~ ), some
-     don't (~ls -l~) and for some arguments the argument is optional
-     (~git --log[=<n>]~)
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-<<type>>
-   #+end_src
-
-** Mutable attributes
-
-   All the options are stored in the hash table ~dict~ where the key
-   is the long option and the value is a value of type ~description~.
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let header_doc      = ref ""
-let description_doc = ref ""
-let footer_doc      = ref ""
-let anon_args_ref   = ref []
-let specs           = ref []
-let dict = Hashtbl.create 67
-   #+end_src
-
-   #+begin_src ocaml :tangle (eval mli)
-val set_header_doc : string -> unit
-val set_description_doc : string -> unit
-val set_footer_doc : string -> unit
-   #+end_src
-
-   Functions to set the header, footer and main description of the
-   documentation provided by the ~help~ function:
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let set_header_doc s = header_doc := s
-let set_description_doc s = description_doc := s
-let set_footer_doc s = footer_doc := s
-   #+end_src
-
-
-   #+begin_src ocaml :tangle (eval mli)
-val anonymous : long_opt -> optional -> documentation -> description
-   #+end_src
-
-   Function to create an anonymous argument.
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let anonymous name opt doc =
-  { short=' ' ; long=name; opt; doc; arg=Without_arg; }
-   #+end_src
-
-** Text formatting functions                                       :noexport:
-
-   Function to print some text such that it fits on the screen
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let output_text t =
-  Format.printf "@[<v 0>";
-  begin
-    match Str.split (Str.regexp "\n") t with
-    | x :: [] ->
-        Format.printf "@[<hov 0>";
-        Str.split (Str.regexp " ") x
-        |> List.iter (fun y -> Format.printf "@[%s@]@ " y) ;
-        Format.printf "@]"
-    | t ->
-        List.iter (fun x ->
-          Format.printf "@[<hov 0>";
-          Str.split (Str.regexp " ") x
-          |> List.iter (fun y -> Format.printf "@[%s@]@ " y) ;
-          Format.printf "@]@;"
-        ) t
-  end;
-  Format.printf "@]"
-   #+end_src
-
-
-   Function to build the short description of the command-line
-   arguments, such as
-
-   Example:
-   #+begin_example
-my_program -b <string> [-h] [-u <float>] -x <string> [--]
-   #+end_example
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let output_short x =
-  match x.short, x.opt, x.arg  with
-  | ' ', Mandatory, _ -> Format.printf "@[%s@]"   x.long
-  | ' ', Optional , _ -> Format.printf "@[[%s]@]" x.long
-  |  _ , Mandatory, Without_arg -> Format.printf "@[-%c@]"   x.short
-  |  _ , Optional , Without_arg -> Format.printf "@[[-%c]@]" x.short
-  |  _ , Mandatory, With_arg arg -> Format.printf "@[-%c %s@]"   x.short arg
-  |  _ , Optional , With_arg arg -> Format.printf "@[[-%c %s]@]" x.short arg
-  |  _ , Mandatory, With_opt_arg arg -> Format.printf "@[-%c [%s]@]"   x.short arg
-  |  _ , Optional , With_opt_arg arg -> Format.printf "@[[-%c [%s]]@]" x.short arg
-   #+end_src
-
-
-   Function to build the long description of the command-line
-   arguments, such as
-
-   Example:
-   #+begin_example
-  -x  --xyz=<string>              Name of the file containing the nuclear
-                                  coordinates in xyz format
-   #+end_example
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let output_long max_width x =
-  let arg =
-    match x.short, x.arg with
-    | ' ' , _                -> x.long
-    |  _  , Without_arg      -> x.long
-    |  _  , With_arg     arg -> Printf.sprintf "%s=%s" x.long arg
-    |  _  , With_opt_arg arg -> Printf.sprintf "%s[=%s]" x.long arg
-  in
-  let long =
-    let l = String.length arg in
-    arg^(String.make (max_width-l) ' ')
-  in
-  Format.printf "@[<v 0>";
-  begin
-    match x.short with
-    | ' '   -> Format.printf "@[%s       @]" long
-    | short -> Format.printf "@[-%c  --%s @]" short long
-  end;
-  Format.printf "@]";
-  output_text  x.doc
-   #+end_src
-
-** Query functions
-
-   #+begin_src ocaml :tangle (eval mli)
-val get       : long_opt -> string option
-val get_bool  : long_opt -> bool
-val anon_args : unit -> string list
-   #+end_src
-
-   | ~get~       | Returns the argument associated with a long option             |
-   | ~get_bool~  | True if the ~Optional~ argument is present in the command-line |
-   | ~anon_args~ | Returns the list of anonymous arguments                        |
-
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let anon_args () = !anon_args_ref
-
-let help () =
-
-  (* Print the header *)
-  output_text !header_doc;
-  Format.printf "@.@.";
-
-  (* Find the anonymous arguments *)
-  let anon =
-    List.filter (fun x -> x.short = ' ') !specs
-  in
-
-  (* Find the options *)
-  let options =
-    List.filter (fun x -> x.short <> ' ') !specs
-    |> List.sort (fun x y -> Char.compare x.short y.short)
-  in
-
-  (* Find column lengths *)
-  let max_width =
-    List.map (fun x ->
-      ( match x.arg with
-        | Without_arg      -> String.length x.long
-        | With_arg arg     -> String.length x.long + String.length arg
-        | With_opt_arg arg -> String.length x.long + String.length arg + 2
-      )
-      + ( if x.opt = Optional then 2 else 0)
-    ) !specs
-    |> List.fold_left max 0
-  in
-
-
-  (* Print usage *)
-  Format.printf "@[<v>@[<v 2>Usage:@,@,@[<hov 4>@[%s@]" Sys.argv.(0);
-  List.iter (fun x -> Format.printf "@ "; output_short x) options;
-  Format.printf "@ @[[--]@]";
-  List.iter (fun x -> Format.printf "@ "; output_short x;) anon;
-  Format.printf "@]@,@]@,";
-
-
-  (* Print arguments and doc *)
-  Format.printf "@[<v 2>Arguments:@,";
-  Format.printf "@[<v 0>" ;
-  List.iter (fun x -> Format.printf "@ "; output_long max_width x) anon;
-  Format.printf "@]@,@]@,";
-
-
-  (* Print options and doc *)
-  Format.printf "@[<v 2>Options:@,";
-
-  Format.printf "@[<v 0>" ;
-  List.iter (fun x -> Format.printf "@ "; output_long max_width x) options;
-  Format.printf "@]@,@]@,";
-
-
-  (* Print footer *)
-  if !description_doc <> "" then
-    begin
-      Format.printf "@[<v 2>Description:@,@,";
-      output_text !description_doc;
-      Format.printf "@,"
-    end;
-
-  (* Print footer *)
-  output_text !footer_doc;
-  Format.printf "@."
-
-
-let get x =
-  try Some (Hashtbl.find dict x)
-  with Not_found -> None
-
-
-let get_bool x = Hashtbl.mem dict x
-   #+end_src
-
-** Specification
-
-   #+begin_src ocaml :tangle (eval mli)
-val set_specs : description list -> unit
-   #+end_src
-
-   Sets the specifications of the current program from a list of
-   ~descrption~ variables.
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let set_specs specs_in =
-  specs := { short = 'h' ;
-             long  = "help" ;
-             doc   = "Prints the help message." ;
-             arg   = Without_arg ;
-             opt   = Optional ;
-           } :: specs_in;
-
-  let cmd_specs =
-    List.filter (fun x -> x.short != ' ') !specs
-    |> List.map (fun { short ; long ; arg ; _ } ->
-      match arg with
-      |  With_arg _ ->
-          (short, long, None, Some (fun x -> Hashtbl.replace dict long x) )
-      |  Without_arg ->
-          (short, long, Some (fun () -> Hashtbl.replace dict long ""), None)
-      |  With_opt_arg _ ->
-          (short, long, Some (fun () -> Hashtbl.replace dict long ""),
-           Some (fun x -> Hashtbl.replace dict long x) )
-    )
-  in
-
-  Getopt.parse_cmdline cmd_specs (fun x -> anon_args_ref := !anon_args_ref @ [x]);
-
-  if (get_bool "help") then
-    (help () ; exit 0);
-
-  (* Check that all mandatory arguments are set *)
-  List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs
-  |> List.iter (fun x ->
-    match get x.long with
-    | Some _ -> ()
-    | None -> failwith ("Error: --"^x.long^" option is missing.")
-  )
-   #+end_src
diff --git a/common/constants.org b/common/constants.org
deleted file mode 100644
index c5aa782..0000000
--- a/common/constants.org
+++ /dev/null
@@ -1,79 +0,0 @@
-#+begin_src elisp tangle: no :results none :exports none
-(setq pwd (file-name-directory buffer-file-name))
-(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
-(setq lib (concat pwd "lib/"))
-(setq testdir (concat pwd "test/"))
-(setq mli (concat lib name ".mli"))
-(setq ml  (concat lib name ".ml"))
-(setq test-ml  (concat testdir name ".ml"))
-(org-babel-tangle)
-#+end_src 
-
-* Constants
-  :PROPERTIES:
-  :header-args: :noweb yes :comments both
-  :END:
-  All constants used in the program.
-  
-** Thresholds
-
-   #+begin_src ocaml :tangle (eval mli)
-val epsilon          : float
-val integrals_cutoff : float
-   #+end_src
-
-   | ~epsilon~          | Value below which a float is considered null. Default is \epsilon = 2.10^{-15} |
-   | ~integrals_cutoff~ | Cutoff value for integrals. Default is \epsilon                                |
-
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let epsilon = 2.e-15
-let integrals_cutoff = epsilon
-   #+end_src
-
-** Mathematical constants 
-
-   #+begin_src ocaml :tangle (eval mli)
-val pi             : float
-val two_pi         : float
-val sq_pi          : float
-val sq_pi_over_two : float
-val pi_inv         : float
-val two_over_sq_pi : float
-   #+end_src
-
-   | ~pi~             | $\pi = 3.141~592~653~589~793~12$ |
-   | ~two_pi~         | $2 \pi$                          |
-   | ~sq_pi~          | $\sqrt{\pi}$                     |
-   | ~sq_pi_over_two~ | $\sqrt{\pi} / 2$                 |
-   | ~pi_inv~         | $1 / \pi$                        |
-   | ~two_over_sq_pi~ | $2 / \sqrt{\pi}$                 |
-    
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let pi             = acos (-1.)
-let two_pi         = 2. *. pi
-let sq_pi          = sqrt pi
-let sq_pi_over_two = sq_pi *. 0.5
-let pi_inv         = 1. /. pi
-let two_over_sq_pi = 2. /. sq_pi
-   #+end_src
-
-** Physical constants
-   
-   #+begin_src ocaml :tangle (eval mli)
-val a0       : float
-val a0_inv   : float
-val ha_to_ev : float
-val ev_to_ha : float
-   #+end_src
-
-   | ~a0~       | Bohr's radius : $a_0 = 0.529~177~210~67(23)$ angstrom |
-   | ~a0_inv~   | $1 / a_0$                                             |
-   | ~ha_to_ev~ | Hartree to eV conversion factor : $27.211~386~02(17)$ |
-   | ~ev_to_ha~ | eV to Hartree conversion factor : 1 / ~ha_to_ev~      |
-    
-   #+begin_src ocaml :tangle (eval ml) :exports none
-let a0 = 0.529_177_210_67
-let a0_inv = 1. /. a0
-let ha_to_ev = 27.211_386_02
-let ev_to_ha = 1. /. ha_to_ev
-   #+end_src
diff --git a/common/lib/charge.ml b/common/lib/charge.ml
index fb5cc13..6ca68ed 100644
--- a/common/lib/charge.ml
+++ b/common/lib/charge.ml
@@ -1,13 +1,5 @@
-
-
-(* This type should be used for all charges in the program (electrons, nuclei,...). *)
-
-
-(* [[file:~/QCaml/common/charge.org::*Type][Type:2]] *)
 type t = float
-(* Type:2 ends here *)
 
-(* [[file:~/QCaml/common/charge.org::*Conversions][Conversions:2]] *)
 external of_float : float -> t = "%identity"
 external to_float : t -> float = "%identity"
 
@@ -23,9 +15,7 @@ let to_string x =
     Printf.sprintf "%f" x
   else 
     "0.0"
-(* Conversions:2 ends here *)
 
-(* [[file:~/QCaml/common/charge.org::*Simple operations][Simple operations:2]] *)
 let gen_op op =
   fun a b ->
   op (to_float a) (to_float b)
@@ -37,9 +27,6 @@ let ( * ) = gen_op ( *. )
 let ( / ) = gen_op ( /. )
 
 let is_null t = t == 0.
-(* Simple operations:2 ends here *)
 
-(* [[file:~/QCaml/common/charge.org::*Printers][Printers:2]] *)
 let pp ppf x = 
   Format.fprintf ppf "@[%s@]" (to_string x)
-(* Printers:2 ends here *)
diff --git a/common/lib/charge.mli b/common/lib/charge.mli
index 7b316a8..2339b36 100644
--- a/common/lib/charge.mli
+++ b/common/lib/charge.mli
@@ -1,15 +1,9 @@
-(* Type
- * 
- *    <<<~Charge.t~>>> *)
+(** This type should be used for all charges in the program (electrons, nuclei,...). *)
 
-(* [[file:~/QCaml/common/charge.org::*Type][Type:1]] *)
 type t
-(* Type:1 ends here *)
 
-(* Conversions *)
+(** Conversions *)
 
-
-(* [[file:~/QCaml/common/charge.org::*Conversions][Conversions:1]] *)
 val of_float : float -> t
 val to_float : t -> float
 
@@ -18,22 +12,16 @@ val to_int   : t -> int
 
 val of_string: string -> t
 val to_string: t -> string
-(* Conversions:1 ends here *)
 
-(* Simple operations *)
+(** Simple operations *)
 
-
-(* [[file:~/QCaml/common/charge.org::*Simple operations][Simple operations:1]] *)
 val ( + ) : t -> t -> t
 val ( - ) : t -> t -> t
 val ( * ) : t -> float -> t
 val ( / ) : t -> float -> t 
 val is_null : t -> bool
-(* Simple operations:1 ends here *)
-
-(* Printers *)
 
 
-(* [[file:~/QCaml/common/charge.org::*Printers][Printers:1]] *)
+(** Printers *)
+
 val pp : Format.formatter -> t -> unit
-(* Printers:1 ends here *)
diff --git a/common/lib/command_line.ml b/common/lib/command_line.ml
index f1081a1..ffcd9b6 100644
--- a/common/lib/command_line.ml
+++ b/common/lib/command_line.ml
@@ -1,17 +1,5 @@
+(** Types *)
 
-
-(* - <<<Short option>>>: in the command line, a dash with a single character
- *   (ex: =ls -l=)
- * - <<<Long option>>>: in the command line, two dashes with a word
- *   (ex: =ls --directory=)
- * - Command-line options can be ~Mandatory~ or ~Optional~
- * - Documentation of the option is used in the help function
- * - Some options require an argument (~ls --ignore="*.ml"~ ), some
- *   don't (~ls -l~) and for some arguments the argument is optional
- *   (~git --log[=<n>]~) *)
-
-
-(* [[file:~/QCaml/common/command_line.org::*Types][Types:2]] *)
 type short_opt     = char
 type long_opt      = string
 type optional      = Mandatory | Optional
@@ -27,51 +15,32 @@ type description = {
   doc  : documentation ;
   arg  : argument ;
 }
-(* Types:2 ends here *)
-
-(* Mutable attributes
- * 
- *    All the options are stored in the hash table ~dict~ where the key
- *    is the long option and the value is a value of type ~description~. *)
 
 
-(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:1]] *)
+
+(** Mutable attributes *)
+
+(** All the options are stored in the hash table ~dict~ where the key
+    is the long option and the value is a value of type ~description~. *)
+
 let header_doc      = ref ""
 let description_doc = ref ""
 let footer_doc      = ref ""
 let anon_args_ref   = ref []
 let specs           = ref []
 let dict = Hashtbl.create 67
-(* Mutable attributes:1 ends here *)
 
 
-
-(* Functions to set the header, footer and main description of the
- * documentation provided by the ~help~ function: *)
-
-
-(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:3]] *)
 let set_header_doc s = header_doc := s
 let set_description_doc s = description_doc := s
 let set_footer_doc s = footer_doc := s
-(* Mutable attributes:3 ends here *)
 
-
-
-(* Function to create an anonymous argument. *)
-
-
-(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:5]] *)
 let anonymous name opt doc =
   { short=' ' ; long=name; opt; doc; arg=Without_arg; }
-(* Mutable attributes:5 ends here *)
 
-(* Text formatting functions                                       :noexport:
- * 
- *    Function to print some text such that it fits on the screen *)
+(** Text formatting functions *)
 
-
-(* [[file:~/QCaml/common/command_line.org::*Text formatting functions][Text formatting functions:1]] *)
+(** Function to print some text such that it fits on the screen *)
 let output_text t =
   Format.printf "@[<v 0>";
   begin
@@ -90,21 +59,15 @@ let output_text t =
         ) t
   end;
   Format.printf "@]"
-(* Text formatting functions:1 ends here *)
 
 
 
+(** Function to build the short description of the command-line arguments, such as
 
-(*    Function to build the short description of the command-line
- *    arguments, such as
- * 
- *    Example:
- *    #+begin_example
- * my_program -b <string> [-h] [-u <float>] -x <string> [--]
- *    #+end_example *)
-
-
-(* [[file:~/QCaml/common/command_line.org::*Text formatting functions][Text formatting functions:2]] *)
+   #+begin_example
+my_program -b <string> [-h] [-u <float>] -x <string> [--]
+   #+end_example
+*)
 let output_short x =
   match x.short, x.opt, x.arg  with
   | ' ', Mandatory, _ -> Format.printf "@[%s@]"   x.long
@@ -115,22 +78,15 @@ let output_short x =
   |  _ , Optional , With_arg arg -> Format.printf "@[[-%c %s]@]" x.short arg
   |  _ , Mandatory, With_opt_arg arg -> Format.printf "@[-%c [%s]@]"   x.short arg
   |  _ , Optional , With_opt_arg arg -> Format.printf "@[[-%c [%s]]@]" x.short arg
-(* Text formatting functions:2 ends here *)
 
 
+(**  Function to build the long description of the command-line arguments, such as
 
-
-(*  Function to build the long description of the command-line
- *  arguments, such as
- * 
- *  Example:
- *  #+begin_example
- * -x  --xyz=<string>              Name of the file containing the nuclear
- *                                 coordinates in xyz format
- *  #+end_example *)
-
-
-(* [[file:~/QCaml/common/command_line.org::*Text formatting functions][Text formatting functions:3]] *)
+  #+begin_example
+ -x  --xyz=<string>              Name of the file containing the nuclear
+                                 coordinates in xyz format
+  #+end_example
+*)
 let output_long max_width x =
   let arg =
     match x.short, x.arg with
@@ -151,17 +107,9 @@ let output_long max_width x =
   end;
   Format.printf "@]";
   output_text  x.doc
-(* Text formatting functions:3 ends here *)
 
 
 
-(* | ~get~       | Returns the argument associated with a long option             |
- * | ~get_bool~  | True if the ~Optional~ argument is present in the command-line |
- * | ~anon_args~ | Returns the list of anonymous arguments                        | *)
-
-
-
-(* [[file:~/QCaml/common/command_line.org::*Query functions][Query functions:2]] *)
 let anon_args () = !anon_args_ref
 
 let help () =
@@ -237,15 +185,8 @@ let get x =
 
 
 let get_bool x = Hashtbl.mem dict x
-(* Query functions:2 ends here *)
 
 
-
-(* Sets the specifications of the current program from a list of
- * ~descrption~ variables. *)
-
-
-(* [[file:~/QCaml/common/command_line.org::*Specification][Specification:2]] *)
 let set_specs specs_in =
   specs := { short = 'h' ;
              long  = "help" ;
@@ -280,4 +221,3 @@ let set_specs specs_in =
     | Some _ -> ()
     | None -> failwith ("Error: --"^x.long^" option is missing.")
   )
-(* Specification:2 ends here *)
diff --git a/common/lib/command_line.mli b/common/lib/command_line.mli
index 12e9a0e..f9658a8 100644
--- a/common/lib/command_line.mli
+++ b/common/lib/command_line.mli
@@ -1,15 +1,71 @@
-(* Types
- * 
- *    #+NAME:type *)
+(** This module is a wrapper around the ~Getopt~ library and helps to
+   define command-line arguments.
+
+   Here is an example of how to use this module.
+   First, define the specification:
+
+  #+begin_src ocaml :tangle no
+let open Command_line in
+begin
+  set_header_doc (Sys.argv.(0) ^ " - One-line description");
+  set_description_doc "Long description of the command.";
+  set_specs
+    [ { short='c'; long="check"; opt=Optional;
+        doc="Checks the input data";
+        arg=Without_arg; };
+
+      { short='b' ; long="basis" ; opt=Mandatory;
+        arg=With_arg "<string>";
+        doc="Name of the file containing the basis set"; } ;
+
+      { short='m' ; long="multiplicity" ; opt=Optional;
+        arg=With_arg "<int>";
+        doc="Spin multiplicity (2S+1). Default is singlet"; } ;
+    ]
+end;
+  #+end_src
+
+   Then, define what to do with the arguments:
+  #+begin_src ocaml :tangle no
+let c =
+  Command_line.get_bool "check"
+in
+
+let basis =
+  match Command_line.get "basis" with
+  | Some x -> x
+  | None -> assert false
+in
+
+let multiplicity =
+  match Command_line.get "multiplicity" with
+  | None -> 1
+  | Some n -> int_of_string n
+in
+  #+end_src
+*)
+
+(** Types *)
 
-(* [[file:~/QCaml/common/command_line.org::type][type]] *)
 type short_opt     = char
+(** In the command line, a dash with a single character (ex: =ls -l=) *)
+  
 type long_opt      = string
+(** In the command line, two dashes with a word (ex: =ls --directory=) *)
+
 type optional      = Mandatory | Optional
+(** Command-line options can be ~Mandatory~ or ~Optional~ *)
+
 type documentation = string
+(** Documentation of the option is used in the help function *)
+  
 type argument      = | With_arg of string
                      | Without_arg
                      | With_opt_arg of string
+(** Some options require an argument (~ls --ignore="*.ml"~ ), some
+    don't (~ls -l~) and for some arguments the argument is optional
+    (~git --log[=<n>]~) *)
+
 
 type description = {
   short: short_opt ;
@@ -18,30 +74,34 @@ type description = {
   doc  : documentation ;
   arg  : argument ;
 }
-(* type ends here *)
 
-(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:2]] *)
+(** Functions to set the header, footer and main description of the
+    documentation provided by the ~help~ function: *)
+
 val set_header_doc : string -> unit
 val set_description_doc : string -> unit
 val set_footer_doc : string -> unit
-(* Mutable attributes:2 ends here *)
 
-(* [[file:~/QCaml/common/command_line.org::*Mutable attributes][Mutable attributes:4]] *)
+
 val anonymous : long_opt -> optional -> documentation -> description
-(* Mutable attributes:4 ends here *)
-
-(* Query functions *)
+(** Function to create an anonymous argument. *)
+
+
+(** Query functions *)
 
 
-(* [[file:~/QCaml/common/command_line.org::*Query functions][Query functions:1]] *)
 val get       : long_opt -> string option
+(** Returns the argument associated with a long option *)
+
 val get_bool  : long_opt -> bool
+(**  True if the ~Optional~ argument is present in the command-line *)
+
 val anon_args : unit -> string list
-(* Query functions:1 ends here *)
-
-(* Specification *)
+(** Returns the list of anonymous arguments *)
 
 
-(* [[file:~/QCaml/common/command_line.org::*Specification][Specification:1]] *)
+(** Specification *)
+
 val set_specs : description list -> unit
-(* Specification:1 ends here *)
+(** Sets the specifications of the current program from a list of
+    ~descrption~ variables. *)
diff --git a/common/lib/constants.ml b/common/lib/constants.ml
index db2233c..34dcb20 100644
--- a/common/lib/constants.ml
+++ b/common/lib/constants.ml
@@ -1,44 +1,22 @@
+(** Thresholds *)
 
-
-(* | ~epsilon~          | Value below which a float is considered null. Default is \epsilon = 2.10^{-15} |
- * | ~integrals_cutoff~ | Cutoff value for integrals. Default is \epsilon                                | *)
-
-
-(* [[file:~/QCaml/common/constants.org::*Thresholds][Thresholds:2]] *)
 let epsilon = 2.e-15
 let integrals_cutoff = epsilon
-(* Thresholds:2 ends here *)
 
 
+(** Mathematical constants  *)
 
-(* | ~pi~             | $\pi = 3.141~592~653~589~793~12$ |
- * | ~two_pi~         | $2 \pi$                          |
- * | ~sq_pi~          | $\sqrt{\pi}$                     |
- * | ~sq_pi_over_two~ | $\sqrt{\pi} / 2$                 |
- * | ~pi_inv~         | $1 / \pi$                        |
- * | ~two_over_sq_pi~ | $2 / \sqrt{\pi}$                 | *)
-
-
-(* [[file:~/QCaml/common/constants.org::*Mathematical constants][Mathematical constants:2]] *)
 let pi             = acos (-1.)
 let two_pi         = 2. *. pi
 let sq_pi          = sqrt pi
 let sq_pi_over_two = sq_pi *. 0.5
 let pi_inv         = 1. /. pi
 let two_over_sq_pi = 2. /. sq_pi
-(* Mathematical constants:2 ends here *)
 
 
+(** Physical constants *)
 
-(* | ~a0~       | Bohr's radius : $a_0 = 0.529~177~210~67(23)$ angstrom |
- * | ~a0_inv~   | $1 / a_0$                                             |
- * | ~ha_to_ev~ | Hartree to eV conversion factor : $27.211~386~02(17)$ |
- * | ~ev_to_ha~ | eV to Hartree conversion factor : 1 / ~ha_to_ev~      | *)
-
-
-(* [[file:~/QCaml/common/constants.org::*Physical constants][Physical constants:2]] *)
 let a0 = 0.529_177_210_67
 let a0_inv = 1. /. a0
 let ha_to_ev = 27.211_386_02
 let ev_to_ha = 1. /. ha_to_ev
-(* Physical constants:2 ends here *)
diff --git a/common/lib/constants.mli b/common/lib/constants.mli
index 0c85ba6..525c77b 100644
--- a/common/lib/constants.mli
+++ b/common/lib/constants.mli
@@ -1,29 +1,46 @@
-(* Thresholds *)
+(** Thresholds *)
 
-
-(* [[file:~/QCaml/common/constants.org::*Thresholds][Thresholds:1]] *)
 val epsilon          : float
+(** Value below which a float is considered null. Default is \epsilon = 2.e-15 *)
+
 val integrals_cutoff : float
-(* Thresholds:1 ends here *)
-
-(* Mathematical constants  *)
+(** Cutoff value for integrals. Default is \epsilon *)
 
 
-(* [[file:~/QCaml/common/constants.org::*Mathematical constants][Mathematical constants:1]] *)
+(** Mathematical constants  *)
+
 val pi             : float
+(** $\pi = 3.141~592~653~589~793~12$ *)
+
 val two_pi         : float
+(** $2 \pi$                          *)
+
 val sq_pi          : float
+(** $\sqrt{\pi}$                     *)
+
 val sq_pi_over_two : float
+(** $\sqrt{\pi} / 2$                 *)
+
 val pi_inv         : float
+(** $1 / \pi$                        *)
+
 val two_over_sq_pi : float
-(* Mathematical constants:1 ends here *)
+(** $2 / \sqrt{\pi}$                 *)
 
-(* Physical constants *)
-   
 
-(* [[file:~/QCaml/common/constants.org::*Physical constants][Physical constants:1]] *)
+
+(** Physical constants *)
+
 val a0       : float
+(** Bohr's radius : $a_0 = 0.529~177~210~67(23)$ angstrom *)
+
 val a0_inv   : float
+(** $1 / a_0$                                             *)
+
 val ha_to_ev : float
+(** Hartree to eV conversion factor : $27.211~386~02(17)$ *)
+
 val ev_to_ha : float
-(* Physical constants:1 ends here *)
+(** eV to Hartree conversion factor : 1 / ~ha_to_ev~      *)
+
+
diff --git a/common/lib/smallarray.ml b/common/lib/smallarray.ml
new file mode 100644
index 0000000..542edbc
--- /dev/null
+++ b/common/lib/smallarray.ml
@@ -0,0 +1,103 @@
+type t = { dim1 : int;
+           dim2 : int;
+           dim3 : int;
+           data : float array;
+         }
+
+let dim1 t = t.dim1
+let dim2 t = t.dim2
+let dim3 t = t.dim3
+
+let data t = t.data
+
+let at t i j k =
+  k + t.dim3*(j + t.dim2*i)
+
+let from t l =
+  let i = l / (t.dim3*t.dim2) in
+  let l' = l - i*t.dim3*t.dim2 in
+  let j = l' / t.dim3 in
+  let k = l' - j * t.dim3
+  in (i,j,k)
+
+let init dim1 dim2 dim3 f =
+  let t = { dim1 ; dim2 ; dim3 ; data = Array.create_float (dim1 * dim2 * dim3) } in
+  let l = ref 0 in
+  for i=0 to dim1-1 do
+    for j=0 to dim2-1 do
+      for k=0 to dim3-1 do
+        t.data.(!l) <- f i j k;
+        incr l
+      done
+    done
+  done;
+  t
+
+let make dim1 dim2 dim3 v =
+  { dim1 ; dim2 ; dim3 ;
+    data = Array.make (dim1 * dim2 * dim3) v
+  }
+
+let get t i j k =
+  Array.get t.data (at t i j k)
+
+let unsafe_get t i j k =
+  Array.unsafe_get t.data (at t i j k)
+
+let set t i j k v =
+  Array.set t.data (at t i j k) v
+
+let unsafe_set t i j k v =
+  Array.unsafe_set t.data (at t i j k) v
+
+let iter3 i j f t =
+  assert (0 <= i && i < t.dim1);
+  assert (0 <= j && j < t.dim2);
+  let start = t.dim3*(j+i*t.dim2) in
+  for k=0 to t.dim3-1 do
+    f (Array.unsafe_get t.data (k+start))
+  done
+
+let iteri3 i j f t =
+  assert (0 <= i && i < t.dim1);
+  assert (0 <= j && j < t.dim2);
+  let start = t.dim3*(j+i*t.dim2) in
+  for k=0 to t.dim3-1 do
+    f k (Array.unsafe_get t.data (k+start))
+  done
+
+let to_array t =
+  Array.init t.dim1 (fun i ->
+    Array.init t.dim2 (fun j ->
+      Array.init t.dim3 (fun k ->
+        unsafe_get t i j k)))
+
+let of_array a =
+  init (Array.length a) (Array.length a.(0)) (Array.length a.(0).(0))
+    (fun i j k -> a.(i).(j).(k) )
+
+let sub1 i t =
+  let n = t.dim2 * t.dim3 in
+  let shift = i * n in
+  { dim1 = 1 ; dim2 = t.dim2 ; dim3 = t.dim3 ;
+    data = Array.init n (fun k -> Array.unsafe_get t.data (shift+k)) }
+
+let sum t =
+  Array.fold_left (+.) 0. t.data
+
+let sum_sub1 k t =
+  let n = t.dim2 * t.dim3 in
+  let shift = k * n in
+  let accu = ref 0. in
+  for i=shift to shift+n-1 do
+    accu := !accu +. Array.unsafe_get t.data i
+  done;
+  !accu
+
+let scale f t =
+  { dim1 = t.dim1 ; dim2 = t.dim2 ; dim3 = t.dim3 ;
+    data = Array.map (fun x -> f *. x) t.data }
+
+let ( .@() ) t (i,j,k) = get t i j k
+let ( .@() <- ) t (i,j,k) v = set t i j k v
+
diff --git a/common/lib/smallarray.mli b/common/lib/smallarray.mli
new file mode 100644
index 0000000..1386a12
--- /dev/null
+++ b/common/lib/smallarray.mli
@@ -0,0 +1,24 @@
+type t = { dim1 : int; dim2 : int; dim3 : int; data : float array; }
+val dim1 : t -> int
+val dim2 : t -> int
+val dim3 : t -> int
+val data : t -> float array
+val at : t -> int -> int -> int -> int
+val from : t -> int -> int * int * int
+val init : int -> int -> int -> (int -> int -> int -> float) -> t
+val make : int -> int -> int -> float -> t
+val get : t -> int -> int -> int -> float
+val unsafe_get : t -> int -> int -> int -> float
+val set : t -> int -> int -> int -> float -> unit
+val unsafe_set : t -> int -> int -> int -> float -> unit
+val iter3 : int -> int -> (float -> unit) -> t -> unit
+val iteri3 : int -> int -> (int -> float -> unit) -> t -> unit
+val to_array : t -> float array array array
+val of_array : float array array array -> t
+val sub1 : int -> t -> t
+val sum : t -> float
+val sum_sub1 : int -> t -> float
+val scale : float -> t -> t
+
+val ( .@() ) : t -> int * int * int -> float
+val ( .@()<- ) : t -> int * int * int -> float -> unit
diff --git a/gaussian_integrals/lib/two_electron_rr_vectorized.ml b/gaussian_integrals/lib/two_electron_rr_vectorized.ml
index 216e817..6f8c4b5 100644
--- a/gaussian_integrals/lib/two_electron_rr_vectorized.ml
+++ b/gaussian_integrals/lib/two_electron_rr_vectorized.ml
@@ -10,6 +10,7 @@ module Cspc = Contracted_shell_pair_couple
 module Po   = Powers
 module Psp  = Primitive_shell_pair
 module Ps   = Primitive_shell
+module Sa   = Smallarray
 module Zp   = Zero_m_parameters
 
 exception NullQuartet
@@ -35,7 +36,7 @@ type four_idx_intermediate =
   center_pq  : Co.axis -> float array array;
   center_pa  : Co.axis -> float array;
   center_qc  : Co.axis -> float array;
-  zero_m_array : float array array array;
+  zero_m_array : Sa.t;
 }
 
 
@@ -54,7 +55,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
 
   let zero_m_array = abcd.zero_m_array in
 
-  let maxm = Array.length zero_m_array - 1 in
+  let maxm = Sa.dim1 zero_m_array - 1 in
 
   let get_xyz angMom = 
     match angMom with
@@ -77,62 +78,48 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
           let xyz = get_xyz angMom_a in
           let am  = Po.decr xyz angMom_a in
           let cab = Co.get xyz center_ab in
-          let result = Array.init (maxm+1-angMom_a.Po.tot) (fun _ -> Array.make_matrix np nq 0.) in
+          let result = Sa.make (maxm+1-angMom_a.Po.tot) np nq 0. in
           let v_am= vrr0_v am in
 
           begin
             if abs_float cab >= cutoff then
-                let expo_b = abcd.expo_b  in
-                Array.iteri (fun m result_m -> 
-                    let v0 = v_am.(m) in
-                    Array.iteri (fun l result_ml ->
-                        let f0 = -. expo_b.(l) *. expo_p_inv.(l) *. cab
-                        and v0_l = v0.(l) 
-                        in
-                        Array.iteri (fun k v0_lk ->
-                            result_ml.(k) <- v0_lk *. f0) v0_l
-                    ) result_m
-                ) result
+              let expo_b = abcd.expo_b  in
+              for m=0 to Sa.dim1 result - 1 do
+                for l=0 to Sa.dim2 result - 1 do
+                  let f0 = -. expo_b.(l) *. expo_p_inv.(l) *. cab in
+                  Sa.iteri3 m l (fun k _ ->
+                    Sa.(result.@(m,l,k) <- v_am.@(m,l,k) *. f0)) result
+                done
+              done
           end;
           let amxyz = Po.get  xyz am in
           if amxyz < 1 then
-              Array.iteri (fun l expo_inv_p_l ->
-                  let center_pq_xyz_l = (center_pq xyz).(l) in
-                  Array.iteri (fun m result_m -> 
-                        let result_ml = result_m.(l) in
-                        let p0  = v_am.(m+1) in
-                        let p0_l = p0.(l)
-                        in
-                        Array.iteri (fun k p0_lk ->
-                            result_ml.(k) <- result_ml.(k)
-                                            +. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk
-                          ) p0_l
-                      ) result
-                ) expo_p_inv
+            Array.iteri (fun l expo_inv_p_l ->
+              let center_pq_xyz_l = (center_pq xyz).(l) in
+              for m=0 to Sa.dim1 result - 1 do
+                Sa.iteri3 m l (fun k result_mlk ->
+                  Sa.(result.@(m,l,k) <- result_mlk
+                                         +. expo_inv_p_l *. center_pq_xyz_l.(k) *. v_am.@(m+1,l,k))
+                ) result
+              done
+            ) expo_p_inv
           else
-              begin
-                let amm = Po.decr xyz am in
-                let amxyz = Util.float_of_int_fast amxyz  in
-                let v_amm = vrr0_v amm in
-                Array.iteri (fun l expo_inv_p_l ->
-                    let f = amxyz *. expo_p_inv.(l) *. 0.5 
-                    and center_pq_xyz_l = (center_pq xyz).(l)
-                    in
-                    Array.iteri (fun m result_m ->
-                        let v1 = v_amm.(m) in
-                        let v1_l = v1.(l) in
-                        let result_ml = result_m.(l) in
-                        let v2 = v_amm.(m+1)  in
-                        let p0  = v_am.(m+1) in
-                        let v2_l = v2.(l)
-                        in
-                        Array.iteri (fun k p0_lk ->
-                            result_ml.(k) <- result_ml.(k) +.
-                                            expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +.
-                                            f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l) 
-                          ) p0.(l)
-                        ) result
-                ) expo_p_inv
+            begin
+              let amm = Po.decr xyz am in
+              let amxyz = Util.float_of_int_fast amxyz  in
+              let v_amm = vrr0_v amm in
+              Array.iteri (fun l expo_inv_p_l ->
+                let f = amxyz *. expo_p_inv.(l) *. 0.5 
+                and center_pq_xyz_l = (center_pq xyz).(l)
+                in
+                for m=0 to Sa.dim1 result - 1 do
+                  Sa.iteri3 m l (fun k result_mlk ->
+                    Sa.(result.@(m,l,k) <- result_mlk +.
+                                     expo_inv_p_l *. center_pq_xyz_l.(k) *. v_am.@(m+1,l,k) +.
+                                     f *. (v_amm.@(m,l,k) +. v_amm.@(m+1,l,k) *. expo_inv_p_l) )
+                  ) result
+                done
+              ) expo_p_inv
               end;
           result 
         in
@@ -142,7 +129,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
   and vrr_v m angMom_a angMom_c =
 
     match (angMom_a.Po.tot, angMom_c.Po.tot) with
-    | (_,0) -> Some (vrr0_v angMom_a).(m)
+    | (_,0) -> Some (vrr0_v angMom_a |> Sa.sub1 m)
     | (_,_) ->
 
       let key =  Zkey.of_powers_six angMom_a angMom_c in
@@ -171,38 +158,36 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
                 match vrr_v m angMom_a cm with
                 | None -> None
                 | Some v1 -> 
+                    Some (Sa.init 1 np nq (fun _ l k ->
+                      Sa.(v1.@(0,l,k)) *. f1.(k) ))
+                      (*
                   begin
                     Some (Array.init np (fun l ->
                       let v1_l = v1.(l) in
                       Array.mapi (fun k f1k -> v1_l.(k) *. f1k) f1
                     ) )
                   end
+                         *)
               else None
             in
 
             let v2 =
               let f2 =
-                Array.init np (fun l ->
+                Sa.init 1 np nq (fun _ l k ->
                     let cpq_l = (center_pq xyz).(l) in
-                    Array.init nq (fun k ->
-                        let x = expo_q_inv.(k) *. cpq_l.(k) in
-                        if (!do_compute) then x
-                        else (if abs_float x > cutoff then do_compute := true ; x)
-                      ) )
+                    let x = expo_q_inv.(k) *. cpq_l.(k) in
+                    if (!do_compute) then x
+                    else (if abs_float x > cutoff then do_compute := true ; x)
+                )
               in
               if (!do_compute) then
                 match vrr_v (m+1) angMom_a cm with
                 | None -> None
                 | Some v2 -> 
                   begin
-                    for l=0 to np-1 do
-                      let f2_l = f2.(l)
-                      and v2_l = v2.(l)
-                      in
-                      for k=0 to nq-1 do
-                        f2_l.(k) <- -. v2_l.(k) *. f2_l.(k)
-                      done
-                    done;
+                    let data_f2 = Sa.data f2
+                    and data_v2 = Sa.data v2 in
+                    Array.iteri (fun i x -> data_f2.(i) <- -. x *. data_f2.(i)) data_v2;
                     Some f2
                   end
               else None
@@ -215,14 +200,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
               | Some v1, None -> Some v1
               | Some v1, Some v2 ->
                 begin
-                  for l=0 to np-1 do
-                    let v1_l = v1.(l)
-                    and v2_l = v2.(l)
-                    in
-                    for k=0 to nq-1 do
-                      v2_l.(k) <- v2_l.(k) +. v1_l.(k)
-                    done
-                  done;
+                  let v1_data = Sa.data v1
+                  and v2_data = Sa.data v2 in
+                  Array.iteri (fun i x -> v2_data.(i) <- v2_data.(i) +. x) v1_data;
                   Some v2
                 end
             in
@@ -245,14 +225,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
                     | None -> None
                     | Some v1 -> 
                       begin
-                        let result = Array.make_matrix np nq 0. in
+                        let result = Sa.make 1 np nq 0. in
                         for l=0 to np-1 do
-                          let v1_l = v1.(l)
-                          and result_l = result.(l)
-                          in
-                          for k=0 to nq-1 do
-                            result_l.(k) <- v1_l.(k) *. f1.(k)
-                          done;
+                           Sa.iteri3 0 l (fun k v1_0lk ->
+                              Sa.(result.@(0,l,k) <- v1_0lk *. f1.(k))) v1
                         done;
                         Some result
                       end
@@ -272,14 +248,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
                     | None -> None
                     | Some v3 -> 
                       begin
-                        let result = Array.make_matrix np nq 0. in
+                        let result = Sa.make 1 np nq 0. in
                         for l=0 to np-1 do
-                          let v3_l = v3.(l)
-                          and result_l = result.(l)
-                          in
-                          for k=0 to nq-1 do
-                            result_l.(k) <- v3_l.(k) *. f2.(k)
-                          done
+                           Sa.iteri3 0 l (fun k v1_0lk ->
+                              Sa.(result.@(0,l,k) <- v1_0lk *. f2.(k))) v3
                         done;
                         Some result
                       end
@@ -292,51 +264,40 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
                 | None, None, Some v3 -> Some v3
                 | Some p1, Some v1, Some v3 ->
                   begin
-                    for l=0 to np-1 do
-                      let v3_l = v3.(l)
-                      and v1_l = v1.(l)
-                      and p1_l = p1.(l)
-                      in
-                      for k=0 to nq-1 do
-                        v3_l.(k) <- p1_l.(k) +. v1_l.(k) +. v3_l.(k)
-                      done
-                    done;
+                    let v3_data = Sa.data v3
+                    and p1_data = Sa.data p1
+                    and v1_data = Sa.data v1 in
+                    Array.iteri (fun i p ->
+                      v3_data.(i) <- p +. v1_data.(i) +. v3_data.(i) ) p1_data;
                     Some v3
                   end
                 | Some p1, Some v1, None ->
                   begin
-                    for l=0 to np-1 do
-                      let v1_l = v1.(l)
-                      and p1_l = p1.(l)
-                      in
-                      for k=0 to nq-1 do
-                        p1_l.(k) <- v1_l.(k) +. p1_l.(k)
-                      done
-                    done;
+                    let p1_data = Sa.data p1
+                    and v1_data = Sa.data v1 in
+                    Array.iteri (fun i p ->
+                      p1_data.(i) <- v1_data.(i) +. p ) p1_data;
                     Some p1
                   end
                 | Some p1, None, Some v3 ->
                   begin
                     for l=0 to np-1 do
-                      let v3_l = v3.(l)
-                      and p1_l = p1.(l)
-                      in
                       for k=0 to nq-1 do
-                        p1_l.(k) <- p1_l.(k) +. v3_l.(k)
+                        Sa.(p1.@(0,l,k) <- p1.@(0,l,k) +. v3.@(0,l,k))
                       done
                     done;
+                    let v3_data = Sa.data v3
+                    and p1_data = Sa.data p1 in
+                    Array.iteri (fun i v ->
+                      p1_data.(i) <- v +. p1_data.(i) ) v3_data;
                     Some p1
                   end
                 | None   , Some v1, Some v3 ->
                   begin
-                    for l=0 to np-1 do
-                      let v3_l = v3.(l)
-                      and v1_l = v1.(l)
-                      in
-                      for k=0 to nq-1 do
-                        v3_l.(k) <- v1_l.(k) +. v3_l.(k)
-                      done
-                    done;
+                    let v3_data = Sa.data v3
+                    and v1_data = Sa.data v1 in
+                    Array.iteri (fun i v ->
+                      v3_data.(i) <- v +. v3_data.(i) ) v1_data;
                     Some v3
                   end
             in
@@ -352,16 +313,13 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
                 begin
                   let p2 =
                     match p2 with
-                    | None -> Array.make_matrix np nq 0.
+                    | None -> Sa.make 1 np nq 0.
                     | Some p2 -> p2
                   in
                   for l=0 to np-1 do
                     let fa = (Util.float_of_int_fast axyz) *. expo_p_inv.(l) *. 0.5 in
-                    let p2_l = p2.(l)
-                    and v_l = v.(l)
-                    in
                     for k=0 to nq-1 do
-                      p2_l.(k) <- p2_l.(k) -. fa *. expo_q_inv.(k) *. v_l.(k)
+                      Sa.(p2.@(0,l,k) <- p2.@(0,l,k) -. fa *. expo_q_inv.(k) *. v.@(0,l,k));
                     done
                   done;
                   Some p2
@@ -498,10 +456,6 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
           *)
   in
 
-  let sum matrix =
-    Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix
-  in
-
   let vrr_v a c = 
     let v = 
       (*
@@ -512,7 +466,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
       vrr_v 0 a c
     in
     match v with 
-      | Some matrix -> sum matrix
+      | Some matrix -> Sa.sum matrix
       | None -> 0.
   in
 
@@ -524,7 +478,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d)
     | 0 ->
       begin
         match (angMom_a.Po.tot, angMom_c.Po.tot) with
-        | (0,0) -> sum zero_m_array.(0)
+        | (0,0) -> Sa.sum zero_m_array
         | (_,_) -> vrr_v angMom_a angMom_c
       end
     | 1 ->
@@ -750,10 +704,7 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell
               | Co.Z -> result.(2)
           in
           let zero_m_array =
-            let result = 
-              Array.init (maxm+1) (fun _ ->
-                  Array.init np (fun _ -> Array.make nq 0. ) )
-            in
+            let result = Sa.make (maxm+1) np nq 0. in
             let empty = Array.make (maxm+1) 0. in
             let center_qc_tmp = Array.init nq (fun cd ->
                 Coordinate.make { Coordinate.
@@ -795,10 +746,8 @@ let contracted_class_shell_pairs ?operator ~zero_m ?schwartz_p ?schwartz_q shell
                 (* Transpose result *)
                 let coef_ab = coef.(ab) in
                 for m=0 to maxm do
-                  let result_m_ab = result.(m).(ab)
-                  in
                   for cd=0 to nq-1 do
-                    result_m_ab.(cd) <- zero_m_array_tmp.(cd).(m) *. coef_ab.(cd)
+                    Sa.(result.@(m,ab,cd) <- zero_m_array_tmp.(cd).(m) *. coef_ab.(cd))
                   done
                 done
               ) sp;
diff --git a/top/lib/install_printers.ml b/top/lib/install_printers.ml
index bf38de5..13f3af2 100644
--- a/top/lib/install_printers.ml
+++ b/top/lib/install_printers.ml
@@ -1,7 +1,6 @@
 (* [[file:~/QCaml/top/install_printers.org::*Intall printers][Intall printers:3]] *)
 let printers =
   [
-    "Common.Charge.pp" ;
     "Common.Coordinate.pp" ;
     "Common.Powers.pp" ;
     "Common.Range.pp" ;