Added Anderson's algorithm

CI/Spindeterminant.ml

@@ -96,6 +96,9 @@ let negate_phase = function
   | None -> None
 
 
+let of_bitstring ?(phase=Phase.Pos) bitstring =
+  Some { bitstring ; phase }
+
 let of_list l = 
   List.rev l
   |> List.fold_left (fun accu p -> creation p accu) vac

CI/Spindeterminant.mli

@@ -26,6 +26,7 @@ val is_none : t -> bool
 val negate_phase : t -> t
 (** Returns a spin-determinant with the phase reversed. *)
 
+
 (** {1 Second quantization operators} *)
 
 val vac : t
@@ -53,6 +54,11 @@ val particles_of : t -> t -> int list
 (** Returns the list of particles in the excitation from one determinant to another. *)
 
 
+(** {1 Creation} *)
+
+val of_bitstring : ?phase:Phase.t -> Z.t -> t
+(** Creates from a bitstring and an optional phase.*)
+
 val of_list : int list -> t
 (** Builds a spin-determinant from a list of orbital indices. If the creation of the
     spin-determinant is not possible because of Pauli's exclusion principle, a [None]

CI/Spindeterminant_space.ml

@@ -0,0 +1,26 @@
+type t =
+{
+  elec_num : int;
+  mo_basis : MOBasis.t;
+  spindets : Spindeterminant.t array;
+}
+
+let fci_of_mo_basis mo_basis elec_num =
+  let mo_num = MOBasis.size mo_basis in
+  let spindets = 
+    Util.bit_permtutations elec_num mo_num 
+    |> List.map (fun b -> Spindeterminant.of_bitstring b)
+    |> Array.of_list
+  in
+  { elec_num ; mo_basis ; spindets }
+
+
+let pp_spindet_space ppf t =
+  Format.fprintf ppf "@[<v>[";
+  Array.iteri (fun i d -> Format.fprintf ppf "@[<h>@[%d@]@;@[%a@]@]@," i Spindeterminant.pp_spindet d) t.spindets;
+  Format.fprintf ppf "]@]"
+
+
+
+
+

MOBasis/MOBasis.ml

@@ -32,6 +32,9 @@ type t =
 }
 
 
+let size t =
+  Mat.dim2 t.mo_coef 
+  
 
 let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = 
   xt_o_x ~x:mo_coef ~o:ao_matrix
@@ -204,3 +207,4 @@ let of_hartree_fock ~frozen_core = function
 | HF.RHF hf -> of_rhf ~frozen_core hf
 | _ -> assert false
 
+

MOBasis/MOBasis.mli

@@ -28,7 +28,6 @@ type t = private
 }
 
 
-
 val make :  ao_basis:AOBasis.t ->
             mo_type:mo_type ->
             mo_class:mo_class array ->
@@ -41,3 +40,6 @@ val of_hartree_fock : frozen_core:bool -> HartreeFock_type.t -> t
 (** Build MOs from a Restricted Hartree-Fock calculation. *)
 
 
+val size : t -> int
+(** Number of molecular orbitals in the basis *)
+

Parallel_mpi/Farm.ml

@@ -5,8 +5,9 @@ let run_sequential f stream =
     try
       let task = Stream.next stream in
       Some (f task) 
-    with Stream.Failure -> None in
-    Stream.from next
+    with Stream.Failure -> None
+  in
+  Stream.from next
 
 
 (* Multi-process functions *)

Utils/Util.ml

@@ -96,6 +96,17 @@ let fact = function
   | i -> fact_memo.(i)
 
 
+let binom n k =
+  (*TODO : slow function *)
+  assert (n > k);
+  let rec aux n k =
+    if k = 0 || k = n then
+      1
+    else
+      aux (n-1) (k-1) + aux (n-1) k
+  in aux n k
+
+
 let rec pow a = function
   | 0 -> 1.
   | 1 -> a
@@ -225,6 +236,21 @@ let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c =
 
 
 
+(** {2 Bitstring functions} *)
+let bit_permtutations m n =
+
+  let rec aux k u rest =
+    if k=1 then
+      List.rev (u :: rest)
+    else
+      let t   = Z.(logor u (u-one)) in
+      let t'  = Z.(t+one) in
+      let t'' = Z.(shift_right ((logand (lognot t) t') - one)) (Z.trailing_zeros u + 1) in
+      aux (k-1) (Z.logor t' t'') (u :: rest)
+  in
+  aux (binom n m) Z.(shift_left one m - one) []
+      
+
 (** {2 Printers} *)
 
 let pp_float_array_size ppf a = 

Utils/Util.mli

@@ -23,6 +23,9 @@ val fact : int -> float
   @raise Invalid_argument for negative arguments or arguments >100.
   *)
 
+val binom : int -> int -> int
+(** Binomial coefficient. [binom n k] {% $= C_n^k$ %}. *)
+
 val pow : float -> int -> float
 (** Fast implementation of the power function for small integer powers *)
 
@@ -103,6 +106,22 @@ val sym_matrix_of_file : string -> Lacaml.D.mat
 (** Reads a symmetric matrix from a file with format "%d %d %f" corresponding to
     [i, j, A.{i,j}]. *)
 
+
+(** {2 Bitstring functions} *)
+
+val bit_permtutations : int -> int -> Z.t list
+(** [bit_permtutations m n] generates the list of all possible [n]-bit
+    strings with [m] bits set to 1.
+    Algorithm adapted from
+    {{:https://graphics.stanford.edu/~seander/bithacks.html#NextBitPermutation}
+    Bit twiddling hacks}.
+    Example:
+{[
+    bit_permtutations 2 4 = [ 0011 ; 0101 ; 0110 ; 1001 ; 1010 ; 1100 ]
+]}
+*)
+
+
 (** {2 Printers} *)
 val pp_float_array_size : Format.formatter -> float array -> unit
 (** Example: