From a9c5fd16a961397b6e3c95dd8da010838b5af4a2 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 24 Jun 2024 17:44:42 +0200
Subject: [PATCH] Added determinant_space

---
 ci/lib/determinant_space.ml   | 343 ++++++++++++++++++++++++++++++++++
 ci/lib/determinant_space.mli  |  74 ++++++++
 linear_algebra/lib/matrix.mli |   2 +-
 top/install_printers.sh       |  47 +++++
 4 files changed, 465 insertions(+), 1 deletion(-)
 create mode 100644 ci/lib/determinant_space.ml
 create mode 100644 ci/lib/determinant_space.mli
 create mode 100755 top/install_printers.sh

diff --git a/ci/lib/determinant_space.ml b/ci/lib/determinant_space.ml
new file mode 100644
index 0000000..94ed591
--- /dev/null
+++ b/ci/lib/determinant_space.ml
@@ -0,0 +1,343 @@
+(** Data structures for storing the determinant space.
+
+If the space is built as the outer product of all {% $\alpha$ %} and {%
+$\beta$ %} determinants, the storage is of type [Spin]. It is sufficient
+to have the arrays of {% $\alpha$ %} and {% $\beta$ %} spindeterminants.
+
+Otherwise, the space is of type [Arbitrary].
+
+*)
+
+open Common
+open Linear_algebra
+
+type arbitrary_space =
+  {
+    det         : int array array ;
+    det_alfa    : Spindeterminant.t array ;
+    det_beta    : Spindeterminant.t array ;
+    index_start : int array;
+  }
+
+
+type determinant_storage =
+  | Arbitrary of arbitrary_space
+  | Spin      of (Spindeterminant.t array * Spindeterminant.t array)
+
+
+type t =
+  {
+    n_alfa       : int ;
+    n_beta       : int ;
+    mo_class     : Mo.Class.t ;
+    mo_basis     : Mo.Basis.t ;
+    determinants : determinant_storage;
+  }
+
+
+module Ss = Spindeterminant_space
+module Electrons = Particles.Electrons
+module Nuclei = Particles.Nuclei
+
+let n_alfa   t = t.n_alfa
+let n_beta   t = t.n_beta
+let mo_class t = t.mo_class
+let mo_basis t = t.mo_basis
+
+
+let size t =
+  match t.determinants with
+  | Spin (a,b)  ->  (Array.length a) * (Array.length b)
+  | Arbitrary a ->
+      let ndet_a = Array.length a.det_alfa in
+      a.index_start.(ndet_a)
+
+
+let determinant_stream t =
+  match t.determinants with
+  | Arbitrary a ->
+    let det_beta = a.det_beta
+    and det_alfa = a.det_alfa
+    and det = a.det in
+    let n_alfa = Array.length det_alfa in
+    let alfa = ref det_alfa.(0)
+    and det_i_alfa = ref det.(0) in
+    let i_alfa = ref 0
+    and k_beta = ref 0
+    in
+    Stream.from (fun _ ->
+        if !i_alfa = n_alfa then None else
+          begin
+            let i_beta = (!det_i_alfa).(!k_beta) in
+            let beta = det_beta.(i_beta) in
+            let result =
+              Some (Determinant.of_spindeterminants (!alfa) beta)
+            in
+            incr k_beta;
+            if !k_beta = Array.length !det_i_alfa then
+              begin
+                k_beta := 0;
+                incr i_alfa;
+                if !i_alfa < n_alfa then
+                  begin
+                    alfa := det_alfa.(!i_alfa);
+                    det_i_alfa := det.(!i_alfa)
+                  end
+              end;
+            result
+          end
+        )
+
+  | Spin (a,b) ->
+    let na = Array.length a
+    and nb = Array.length b in
+    let i = ref 0
+    and j = ref 0 in
+    Stream.from (fun _ ->
+        if !j < nb then
+          let result =
+            Determinant.of_spindeterminants a.(!i) b.(!j)
+          in
+          incr i;
+          if !i = na then (i := 0 ; incr j);
+          Some result
+        else
+          None)
+
+
+let determinants t = t.determinants
+
+
+let determinants_array t =
+  let s = determinant_stream t in
+  Array.init (size t) (fun _ -> Stream.next s)
+
+
+(*
+let determinant t i =
+  let alfa, beta =
+    match t.determinants with
+    | Arbitrary a ->
+      let i_alfa =
+        let index_start = a.index_start in
+        let rec loop i_alfa =
+          if index_start.(i_alfa) <= i then
+            (loop [@tailcall]) (i_alfa+1)
+          else i_alfa
+        in loop 0
+      in
+      let i_beta = i - a.index_start.(i_alfa) in
+      let alfa = a.det_alfa.(i_alfa) in
+      let beta = a.det_beta.(i_beta) in
+      alfa, beta
+
+    | Spin (a,b) ->
+      let nb = Array.length b in
+      let k = i / nb in
+      let j = i - k * nb in
+      a.(j), b.(k)
+
+    in
+    Determinant.of_spindeterminants alfa beta
+
+*)
+
+let fock_diag det_space det =
+
+  let alfa_list =
+    Determinant.alfa det
+    |> Spindeterminant.to_list
+  in
+  let beta_list =
+    Determinant.beta det
+    |> Spindeterminant.to_list
+  in
+  let mo_basis = mo_basis det_space in
+  let mo_num = Mo.Basis.size mo_basis in
+  let one_e_ints  = Mo.Basis.one_e_ints mo_basis
+  and two_e_ints  = Mo.Basis.two_e_ints mo_basis
+  in
+  let two_e i j k l = Four_idx_storage.get_phys two_e_ints i j k l in
+  let build_array list1 list2 =
+    let result = Array.make (mo_num+1) 0. in
+
+    (* Occupied *)
+    List.iter (fun i ->
+        let x = one_e_ints%:(i,i) in
+        result.(i) <- result.(i) +. x;
+        result.(0) <- result.(0) +. x;
+        List.iter (fun j ->
+            if j <> i then
+              begin
+                let x = two_e i j i j -. two_e i j j i in
+                result.(i) <- result.(i) +. x;
+                result.(0) <- result.(0) +. 0.5 *. x;
+              end
+          ) list1;
+        List.iter (fun j ->
+              begin
+                let x = two_e i j i j in
+                result.(i) <- result.(i) +. x;
+                result.(0) <- result.(0) +. 0.5 *. x;
+              end
+          ) list2;
+      ) list1;
+
+    (* Virtuals*)
+    List.iter (fun i ->
+        if result.(i) = 0. then
+          begin
+            let x = one_e_ints%:(i,i) in
+            result.(i) <- result.(i) +. x;
+            List.iter (fun j ->
+                let x = two_e i j i j -. two_e i j j i in
+                result.(i) <- result.(i) +. x;
+              ) list1;
+            List.iter (fun j ->
+                  begin
+                    let x = two_e i j i j in
+                    result.(i) <- result.(i) +. x;
+                  end
+              ) list2;
+          end
+     ) (Util.list_range 1 mo_num);
+     result
+  in
+  let alfa, beta =
+    build_array alfa_list beta_list,
+    build_array beta_list alfa_list
+  in
+  let e = alfa.(0) +. beta.(0) in
+  alfa.(0) <- e;
+  beta.(0) <- e;
+  alfa, beta
+
+
+
+let spin_of_mo_basis mo_basis f =
+
+  let s = Mo.Basis.simulation mo_basis in
+  let e = Simulation.electrons s in
+  let n_alfa = Electrons.n_alfa e
+  and n_beta = Electrons.n_beta e in
+  let det_a = f n_alfa
+  and det_b = f n_beta
+  in
+  let mo_class = Ss.mo_class det_a in
+  let determinants =
+    let det_alfa = Ss.spin_determinants det_a
+    and det_beta = Ss.spin_determinants det_b
+    in
+    let n_det_beta = Array.length det_beta in
+    let n_det_alfa = Array.length det_alfa in
+
+    let ndet = n_det_alfa * n_det_beta in
+    Format.printf "Number of determinants : %d %d %d\n%!"
+        n_det_alfa n_det_beta ndet;
+    Spin (det_alfa, det_beta)
+  in
+  { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants }
+
+
+(*
+let arbitrary_of_mo_basis mo_basis f =
+
+  let s = Mo.Basis.simulation mo_basis in
+  let e = Simulation.electrons s in
+  let n_alfa = Electrons.n_alfa e
+  and n_beta = Electrons.n_beta e in
+  let det_a = f n_alfa
+  and det_b = f n_beta
+  in
+  let mo_class = Ss.mo_class det_a in
+  let determinants =
+    let det_alfa = Ss.spin_determinants det_a
+    and det_beta = Ss.spin_determinants det_b
+    in
+    let n_det_beta = Array.length det_beta in
+    let n_det_alfa = Array.length det_alfa in
+
+    let det = Array.make n_det_alfa
+        (Array.init n_det_beta (fun i -> i))
+    in
+    let index_start = Array.init (n_det_alfa+1) (fun i -> i*n_det_beta) in
+    let ndet = (index_start.(n_det_alfa)) in
+
+    Format.printf "Number of determinants : %d %d %d\n%!"
+        n_det_alfa n_det_beta ndet;
+    Arbitrary {
+      det_alfa ; det_beta ; det ; index_start
+    }
+  in
+  { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants }
+*)
+
+
+
+let cas_of_mo_basis mo_basis ~frozen_core n m =
+  let f n_alfa =
+    Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m
+  in
+  spin_of_mo_basis mo_basis f
+
+
+let fci_of_mo_basis mo_basis ~frozen_core =
+  let f n_alfa =
+    Ss.fci_of_mo_basis mo_basis ~frozen_core n_alfa
+  in
+  spin_of_mo_basis mo_basis f
+
+
+let fci_f12_of_mo_basis mo_basis ~frozen_core mo_num =
+  let s = Mo.Basis.simulation mo_basis in
+  let e = Simulation.electrons s in
+  let n_alfa = Electrons.n_alfa e
+  and n_beta = Electrons.n_beta e in
+  let n_core = Mo.Frozen_core.num_mos frozen_core in
+  let n, m =
+      (n_alfa + n_beta - n_core),
+      (mo_num - n_core)
+  in
+  let f n_alfa =
+    Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m
+  in
+  let r =
+    spin_of_mo_basis mo_basis f
+  in
+  { r with mo_class =
+    Mo.Class.to_list r.mo_class
+    |> List.rev_map (fun i ->
+      match i with
+      | Mo.Class.Virtual i when i > mo_num -> Mo.Class.Auxiliary i
+      | i -> i)
+    |> List.rev
+    |> Mo.Class.of_list }
+
+
+let cas_f12_of_mo_basis mo_basis ~frozen_core n m mo_num =
+  let f n_alfa =
+    Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m
+  in
+  let r =
+    spin_of_mo_basis mo_basis f
+  in
+  { r with mo_class =
+    Mo.Class.to_list r.mo_class
+    |> List.rev_map (fun i ->
+      match i with
+      | Mo.Class.Virtual i when i > mo_num -> Mo.Class.Auxiliary i
+      | i -> i)
+    |> List.rev
+    |> Mo.Class.of_list
+    }
+
+
+
+let pp ppf t =
+  Format.fprintf ppf "@[<v 2>[ ";
+  let i = ref 0 in
+  determinant_stream t
+  |> Stream.iter (fun d -> Format.fprintf ppf "@[<v>@[%8d@]@;@[%a@]@]@;" !i
+                    Determinant.pp  d; incr i) ;
+  Format.fprintf ppf "]@]"
+
diff --git a/ci/lib/determinant_space.mli b/ci/lib/determinant_space.mli
new file mode 100644
index 0000000..69bb11c
--- /dev/null
+++ b/ci/lib/determinant_space.mli
@@ -0,0 +1,74 @@
+(**
+The determinant space in which we solve the Schrodinger equation.
+*)
+
+type t
+
+type arbitrary_space =
+  {
+    det         : int array array ;
+    det_alfa    : Spindeterminant.t array ;
+    det_beta    : Spindeterminant.t array ;
+    index_start : int array;
+  }
+
+
+type determinant_storage =
+  | Arbitrary of arbitrary_space
+  | Spin      of (Spindeterminant.t array * Spindeterminant.t array)
+
+
+(** {1 Accessors} *)
+
+val n_alfa : t -> int
+(** Number of {% $\alpha$ %} electrons in the {% $\alpha$ %} MOs. *)
+
+val n_beta : t -> int
+(** Number of {% $\beta$ %} electrons in the {% $\beta$ %} MOs. *)
+
+val mo_class : t -> Mo.Class.t
+(** The MO classes used to generate the space. *)
+
+val mo_basis : t -> Mo.Basis.t
+(** The MO basis on which the determinants are expanded. *)
+
+val determinants : t -> determinant_storage
+(** All the determinants belonging to the space. *)
+
+val determinants_array : t -> Determinant.t array
+(** All the determinants belonging to the space. *)
+
+val determinant_stream : t -> Determinant.t Stream.t
+(** All the determinants belonging to the space, as a stream. *)
+
+val size : t -> int
+(** Size of the determinant space *)
+
+val fock_diag : t -> Determinant.t -> float array * float array
+(** Returns the diagonal of the {% $\alpha$ %} and {% $\beta$ %} Fock matrices.
+    The zero elements contain the energy of the determinant.
+*)
+
+
+val fci_of_mo_basis : Mo.Basis.t -> frozen_core:Mo.Frozen_core.t -> t
+(** Creates a space of all possible ways to put [n_alfa] electrons in the {% $\alpha$ %}
+    [Active] MOs and [n_beta] electrons in the {% $\beta$ %} [Active] MOs.
+    All other MOs are untouched.
+*)
+
+val cas_of_mo_basis : Mo.Basis.t -> frozen_core:Mo.Frozen_core.t -> int -> int -> t
+(** Creates a CAS(n,m) space of determinants.  *)
+
+val fci_f12_of_mo_basis : Mo.Basis.t -> frozen_core:Mo.Frozen_core.t -> int -> t
+(** Creates the active space to perform a FCI-F12 with an
+    auxiliary basis set. *)
+
+val cas_f12_of_mo_basis : Mo.Basis.t -> frozen_core:Mo.Frozen_core.t -> int -> int -> int -> t
+(** [cas_of_mo_basis mo_basis m n mo_num] Creates a CAS(n,m) space
+  of determinants with an auxiliary basis set defined as the MOs from
+  [mo_num+1] to [Mo.Basis.size mo_basis].
+*)
+
+(** {2 Printing} *)
+
+val pp : Format.formatter -> t -> unit
diff --git a/linear_algebra/lib/matrix.mli b/linear_algebra/lib/matrix.mli
index 3e443fe..d2f8e32 100644
--- a/linear_algebra/lib/matrix.mli
+++ b/linear_algebra/lib/matrix.mli
@@ -337,7 +337,7 @@ val sysv : b:('a,'b) t -> ('a,'a) t -> ('a,'b) t
 (** Solves %{ $AX=B$ %} when A is symmetric *)
   
 val (%:) : ('a,'b) t -> (int*int) -> float
-(** [t%.(i,j)] returns the element at i,j. *)
+(** [t%:(i,j)] returns the element at i,j. *)
 
 val set : ('a,'b) t -> int -> int -> float -> unit
 (** [set t i j v] sets the (i,j)-th element to v *)
diff --git a/top/install_printers.sh b/top/install_printers.sh
new file mode 100755
index 0000000..d921540
--- /dev/null
+++ b/top/install_printers.sh
@@ -0,0 +1,47 @@
+#!/bin/bash
+
+function print_modules ()
+{
+    MODULES=$(grep "val pp " ../*/lib/*.mli \
+                  | cut -d ":" -f 1 \
+                  | sed 's|^../||' \
+                  | sed 's|\.mli||' \
+                  | sed 's|/lib/| |' \
+                  | sed 's/[^ ]*/\u&/g' \
+                  | sed 's| |.|' \
+                  | grep -v Top.Install_printers )
+    
+    for MODULE in $MODULES ; do
+        M=($(echo $MODULE | sed 's|\.| |'))
+        if [[ "${M[0]}" == "${M[1]}" ]] ; then
+            echo "    \"${M[1]}.pp\" ;"
+        else
+            echo "    \"${M[0]}.${M[1]}.pp\" ;"
+        fi
+    done
+}    
+
+cat << EOF
+let printers =
+  [
+$(print_modules)
+  ]
+
+let eval_exn str =
+  let lexbuf = Lexing.from_string str in
+  let phrase = !Toploop.parse_toplevel_phrase lexbuf in
+  Toploop.execute_phrase false Format.err_formatter phrase
+
+
+let rec install_printers = function
+  | [] -> eval_exn "#require \"lacaml.top\";;"
+  | printer :: printers ->
+      let cmd = Printf.sprintf "#install_printer %s;;" printer in
+      eval_exn cmd && install_printers printers
+
+let () =
+  if not (install_printers printers) then
+    Format.eprintf "Problem installing QCaml-printers@."
+
+EOF
+