Added CI

ci/lib/determinant.ml

@@ -0,0 +1,116 @@
+open Common
+
+type t =
+{
+  alfa : Spindeterminant.t ;
+  beta : Spindeterminant.t ;
+}
+
+type hole = int
+type particle = int
+
+
+
+
+let alfa t = t.alfa
+
+let beta  t = t.beta
+
+let vac n =
+  {
+    alfa = Spindeterminant.vac n;
+    beta = Spindeterminant.vac n;
+  }
+
+let phase t =
+  match Spindeterminant.(phase t.alfa, phase t.beta) with
+  | Phase.Pos, Phase.Pos
+  | Phase.Neg, Phase.Neg  -> Phase.Pos
+  | _ -> Phase.Neg
+
+
+let of_spindeterminants a b =
+  {
+    alfa = a ;
+    beta  = b
+  }
+
+
+let is_none t = Spindeterminant.(is_none t.alfa || is_none t.beta)
+
+
+let negate_phase t =
+  { t with alfa = Spindeterminant.negate_phase t.alfa }
+
+let set_phase p t =
+  { alfa = Spindeterminant.set_phase p t.alfa ;
+    beta = Spindeterminant.set_phase Phase.Pos t.beta
+  }
+
+
+let excitation_level_alfa t t' =
+  Spindeterminant.excitation_level t.alfa t'.alfa
+
+let excitation_level_beta t t' =
+  Spindeterminant.excitation_level t.beta t'.beta
+
+let excitation_levels t t' =
+  excitation_level_alfa t t', excitation_level_beta t t'
+
+
+let excitation_level t t' =
+  (excitation_level_alfa t t') + (excitation_level_beta t t')
+
+
+let of_lists n a b =
+  let alfa = Spindeterminant.of_list n a
+  and beta  = Spindeterminant.of_list n b
+  in of_spindeterminants alfa beta
+
+
+let to_lists t =
+  Spindeterminant.to_list t.alfa,
+  Spindeterminant.to_list t.beta
+
+
+let creation spin p t =
+  match spin with
+  | Spin.Alfa -> { t with alfa = Spindeterminant.creation p t.alfa }
+  | Spin.Beta -> { t with beta = Spindeterminant.creation p t.beta }
+
+
+let annihilation spin h t =
+  match spin with
+  | Spin.Alfa -> { t with alfa = Spindeterminant.annihilation h t.alfa }
+  | Spin.Beta -> { t with beta = Spindeterminant.annihilation h t.beta }
+
+
+let single_excitation spin h p t =
+  assert (h <> p);
+  match spin with
+  | Spin.Alfa -> { t with alfa = Spindeterminant.single_excitation h p t.alfa }
+  | Spin.Beta -> { t with beta = Spindeterminant.single_excitation h p t.beta }
+
+
+let double_excitation spin h p spin' h' p' t =
+  assert (h <> p);
+  assert (h' <> p');
+  match spin, spin' with
+  | Spin.(Alfa, Beta)  ->        { alfa = Spindeterminant.single_excitation h  p      t.alfa ;
+                                   beta = Spindeterminant.single_excitation h' p'     t.beta }
+  | Spin.(Beta, Alfa)  ->        { beta = Spindeterminant.single_excitation h  p      t.beta ;
+                                   alfa = Spindeterminant.single_excitation h' p'     t.alfa }
+  | Spin.(Alfa, Alfa)  -> { t with alfa = Spindeterminant.double_excitation h p h' p' t.alfa }
+  | Spin.(Beta, Beta)  -> { t with beta = Spindeterminant.double_excitation h p h' p' t.beta }
+
+
+let compare = compare
+
+
+let pp ppf t =
+  Format.fprintf ppf "@[<v>@[phase:%a@]@;@[a:%a@]@;@[b:%a@]@]@."
+    Phase.pp (phase t)
+    Spindeterminant.pp t.alfa
+    Spindeterminant.pp t.beta
+
+

ci/lib/determinant.mli

@@ -0,0 +1,90 @@
+(** A Slater determinant is expressed as a Waller-Hartree double determinant:
+{% $$
+D(\mathbf{R}) = D_\alpha(\mathbf{R_\alpha}) \times D_\beta(\mathbf{R_\beta}) 
+$$ %}
+The {% $\alpha$ %} and {% $\beta$ %} determinants are of type [Spindeterminant.t].
+*)
+
+open Common
+
+type t
+type hole = int
+type particle = int
+
+(** {1 Accessors} *)
+
+val alfa : t -> Spindeterminant.t
+(** Get the {% $\alpha$ %} spin-determinant. *)
+
+val beta : t -> Spindeterminant.t
+(** Get the {% $\beta$ %} spin-determinant. *)
+
+val phase : t -> Phase.t
+(** Get the phase of the Slater determinant, the product of the phases of the
+    spin-determinants.
+*)
+
+val is_none : t -> bool
+(** Tests if a Determinant is [None]. *)
+
+
+(** {1 Second quantization operators} *)
+
+val vac : int -> t
+(** Vacuum state, [vac = Some ]{% $|\rangle$ %}. The integer parameter is the size of the
+    MO basis set. *)
+
+val creation : Spin.t -> particle -> t -> t
+(** [creation spin p] is the creation operator {% $a^\dagger_p$ %}. *)
+
+val annihilation : Spin.t -> hole -> t -> t
+(** [annihilation spin h] is the annihilation operator {% $a_h$ %}. *)
+
+val single_excitation : Spin.t -> hole -> particle -> t -> t
+(** Single excitation operator {% $T_h^p = a^\dagger_p a_h$ %}. *)
+
+val double_excitation : Spin.t -> hole -> particle -> Spin.t -> hole -> particle -> t -> t
+(** Double excitation operator {% $T_{hh'}^{pp'} = a^\dagger_p a^\dagger_{p'} a_{h'} a_h$ %}. *)
+
+val excitation_level_alfa : t -> t -> int
+(** Returns the excitation level between two determinants in the {% $\alpha$ %} spin. *)
+
+val excitation_level_beta : t -> t -> int
+(** Returns the excitation level between two determinants in the {% $\beta$ %} spin. *)
+
+val excitation_levels : t -> t -> int*int
+(** Returns levels of excitation between two determinants in {% $\alpha$ %} and
+    {% $\beta$ %} spins as a pair. *)
+
+val excitation_level : t -> t -> int
+(** Returns excitation level between two determinants. *)
+
+val to_lists : t -> int list * int list
+(** Converts a Slater determinant to a pair of lists of orbital indices.  *)
+
+
+(** {1 Creators} *)
+
+val of_spindeterminants : Spindeterminant.t -> Spindeterminant.t -> t
+(** Creates a Slater determinant from an {% $\alpha$ %} and a {% $\beta$ %}
+    [Spindeterminant.t].
+*)
+
+val of_lists : int -> int list -> int list -> t
+(** Creates a Slater determinant from a two lists of orbital indices.
+    The integer parameter is the size of the MO basis set. *)
+
+val negate_phase : t -> t
+(** Returns the same determinant with the phase negated. *)
+
+val set_phase : Phase.t -> t -> t
+(** Returns the same determinant with the phase set to [p]. *)
+
+val compare : t -> t -> int
+(** Comparison function for sorting *)
+
+
+(** {1 Printers} *)
+
+val pp : Format.formatter -> t -> unit
+

ci/lib/dune

@@ -0,0 +1,13 @@
+(library
+ (name ci)
+ (public_name qcaml.ci)
+ (synopsis "Configuration interaction")
+ (libraries
+  qcaml.common
+  qcaml.particles
+  qcaml.mo
+  )
+  (instrumentation (backend landmarks))
+ 
+)
+

ci/lib/excitation.ml

@@ -0,0 +1,157 @@
+open Common
+
+type single_exc =
+  {
+    hole     : int    ;
+    particle : int    ;
+    spin     : Spin.t ;
+  }
+
+type t =
+| Identity  of Phase.t
+| Single    of Phase.t * single_exc
+| Double    of Phase.t * single_exc * single_exc
+| Triple    of Phase.t * single_exc * single_exc * single_exc
+| Multiple  of Phase.t * single_exc list
+
+
+let single_of_spindet t t' =
+  assert (Spindeterminant.excitation_level t t' = 1);
+  let d  = Spindeterminant.bitstring t
+  and d' = Spindeterminant.bitstring t'
+  in
+  let tmp = Bitstring.logxor d d' in
+  let shift_left_one = Bitstring.(shift_left_one (numbits tmp)) in
+  let hole_z     = Bitstring.logand (Spindeterminant.bitstring t ) tmp
+  and particle_z = Bitstring.logand (Spindeterminant.bitstring t') tmp
+  in
+  let hole     = 1 + Bitstring.trailing_zeros hole_z
+  and particle = 1 + Bitstring.trailing_zeros particle_z
+  in
+  (* Phase calculation *)
+  let low, high =
+    if particle > hole then hole, particle
+    else particle, hole
+  in
+  let mask =
+    let h = high-1 in
+    let l = low in
+    let mask_up = shift_left_one h |> Bitstring.minus_one
+    and mask_dn = Bitstring.plus_one @@ Bitstring.lognot (shift_left_one l)
+    in Bitstring.logand mask_up mask_dn
+  in
+  let phase =
+    Phase.multiply (Phase.multiply (Spindeterminant.phase t) (Spindeterminant.phase t'))
+      (Phase.of_nperm (Bitstring.popcount @@ Bitstring.logand d mask ))
+  in
+  (hole, particle, phase)
+
+
+let single_of_det t t' =
+  assert Determinant.(beta t = beta t' || alfa t = alfa t');
+  if Determinant.(beta t = beta t') then
+    let hole, particle, phase =
+      single_of_spindet (Determinant.alfa t) (Determinant.alfa t')
+    in
+    Single (phase, { hole ; particle ; spin=Spin.Alfa })
+  else
+    let hole, particle, phase =
+      single_of_spindet (Determinant.beta t) (Determinant.beta t')
+    in
+    Single (phase, { hole ; particle ; spin=Spin.Beta })
+
+
+let multiple_of_spindet t t' =
+  let holes     = Spindeterminant.holes_of t t'
+  and particles = Spindeterminant.particles_of t t'
+  in
+  let t'' =
+    List.fold_left (fun accu h -> Spindeterminant.annihilation h accu) t holes
+  in
+  let t'' =
+    List.fold_left (fun accu p -> Spindeterminant.creation p accu) t'' particles
+  in
+  assert (t' = t'' || t' = Spindeterminant.negate_phase t'');
+  let phase =
+    if Spindeterminant.phase t' = Spindeterminant.phase t'' then
+      Phase.Pos
+    else
+      Phase.Neg
+  in
+  (phase, List.rev @@ List.rev_map2 (fun hole particle -> (hole, particle)) holes (List.rev particles) )
+
+
+let double_of_spindet t t' =
+  match multiple_of_spindet t t' with
+  | (phase, (h1,p1)::(h2,p2)::[]) -> (h1, p1, h2, p2, phase)
+  | _ -> invalid_arg "t and t' are not doubly excited"
+
+
+let triple_of_spindet t t' =
+  match multiple_of_spindet t t' with
+  | (phase, (h1,p1)::(h2,p2)::(h3,p3)::[]) -> (h1, p1, h2, p2, h3, p3, phase)
+  | _ -> invalid_arg "t and t' are not doubly excited"
+
+
+let multiple_of_det t t' =
+  let pa, a =
+    multiple_of_spindet (Determinant.alfa t) (Determinant.alfa t')
+  and pb, b =
+    multiple_of_spindet (Determinant.beta t) (Determinant.beta t')
+  in
+  let phase = Phase.multiply pa pb in
+  Multiple (phase, List.concat [
+    List.rev @@ List.rev_map (fun (hole, particle) -> { hole ; particle ; spin=Spin.Alfa }) a ;
+    List.rev @@ List.rev_map (fun (hole, particle) -> { hole ; particle ; spin=Spin.Beta }) b ])
+
+
+let double_of_det t t' =
+  match multiple_of_det t t' with
+  | Multiple (phase, [e1 ; e2]) -> Double (phase, e1, e2)
+  | _ -> assert false
+
+
+let triple_of_det t t' =
+  match multiple_of_det t t' with
+  | Multiple (phase, [e1 ; e2 ; e3]) -> Triple (phase, e1, e2, e3)
+  | _ -> assert false
+
+
+let of_det t t' =
+  match Determinant.excitation_level t t' with
+  | 0 -> if Determinant.phase t = Determinant.phase t' then
+          Identity Phase.Pos
+        else
+          Identity Phase.Neg
+  | 1 -> single_of_det t t'
+  | 2 -> double_of_det t t'
+  | 3 -> triple_of_det t t'
+  | _ -> multiple_of_det t t'
+
+let pp_s_exc ppf t =
+  Format.fprintf ppf "@[T^{%s}_{%d->%d}@]"
+    (match t.spin with
+      | Spin.Alfa -> "\\alpha"
+      | Spin.Beta -> "\\beta " )
+    t.hole t.particle
+
+let pp ppf t =
+  let phase, l =
+    match t with
+    | Identity p -> p, []
+    | Single (p,x) -> p, x::[]
+    | Double (p,x,y) -> p, x::y::[]
+    | Triple (p,x,y,z) -> p, x::y::z::[]
+    | Multiple (p,l) -> p, l
+  in
+  Format.fprintf ppf "@[%c"
+    (match phase with
+      | Phase.Pos -> '+'
+      | Phase.Neg -> '-' );
+  List.iter (fun x -> Format.fprintf ppf "@[T^{%s}_{%d->%d}@]"
+                      (match x.spin with
+                        | Spin.Alfa -> "\\alpha"
+                        | Spin.Beta -> "\\beta " )
+                      x.hole x.particle) l;
+  Format.fprintf ppf "@]"
+

ci/lib/phase.ml

@@ -0,0 +1,31 @@
+type t =
+| Pos
+| Neg
+
+let of_nperm nperm =
+  if (nperm land 1) = 1 then Neg
+  else Pos
+
+let to_nperm = function
+| Pos -> 0
+| Neg -> 1
+
+let multiply t t' =
+  match t, t' with
+  | Pos, Pos
+  | Neg, Neg -> Pos
+  | Pos, Neg
+  | Neg, Pos -> Neg
+
+let neg = function
+  | Pos -> Neg
+  | Neg -> Pos
+
+let add_nperm phase = function
+  | 0 -> phase
+  | nperm -> multiply phase (of_nperm nperm)
+
+let pp ppf = function
+  | Pos -> Format.fprintf ppf "@[<h>+1@]"
+  | Neg -> Format.fprintf ppf "@[<h>-1@]"
+

ci/lib/phase.mli

@@ -0,0 +1,23 @@
+type t =
+| Pos
+| Neg
+
+val of_nperm : int -> t
+(** Returns the phase obtained by a given number of permuations. *)
+
+val to_nperm : t -> int
+(** Converts the phase to [1] or [0] permutations. *)
+
+val multiply : t -> t -> t
+(** Multiply an existing phase by another phase. *)
+
+val add_nperm : t -> int -> t
+(** Add to an existing phase a given number of permutations. *)
+
+val neg : t -> t
+(** Negate the phase. *)
+
+(** {1 Printers} *)
+
+val pp : Format.formatter -> t -> unit
+

ci/lib/spindeterminant.ml

@@ -0,0 +1,144 @@
+open Common
+
+type s =
+  {
+    bitstring : Bitstring.t;
+    phase     : Phase.t ;
+  }
+
+type t = s option
+type hole = int
+type particle = int
+
+let phase = function
+  | Some s -> s.phase
+  | None -> Phase.Pos
+
+
+let is_none = function
+  | None -> true
+  | _ -> false
+
+
+let bitstring = function
+  | Some s -> s.bitstring
+  | None -> invalid_arg "Spindeterminant is None"
+
+
+let vac n =
+  Some  { bitstring = Bitstring.zero n;
+          phase     = Phase.Pos; }
+
+
+let creation p = function
+  | None -> None
+  | Some spindet ->
+    let i = pred p in
+    if Bitstring.testbit spindet.bitstring i then
+      None
+    else
+      begin
+        let numbits = Bitstring.numbits spindet.bitstring in
+        let x = Bitstring.shift_left_one numbits i in
+        let bitstring = Bitstring.logor spindet.bitstring x in
+        let mask = Bitstring.minus_one x in
+        let r = Bitstring.logand bitstring mask in
+        let phase = Phase.add_nperm spindet.phase (Bitstring.popcount r) in
+        Some { bitstring ; phase }
+      end
+
+
+let annihilation h = function
+  | None -> None
+  | Some spindet ->
+    let i = pred h in
+    if not (Bitstring.testbit spindet.bitstring i) then
+      None
+    else
+      begin
+        let numbits = Bitstring.numbits spindet.bitstring in
+        let x = Bitstring.shift_left_one numbits i in
+        let mask = Bitstring.minus_one x in
+        let r = Bitstring.logand spindet.bitstring mask in
+        let phase = Phase.add_nperm spindet.phase (Bitstring.popcount r) in
+        let bitstring = Bitstring.logand spindet.bitstring (Bitstring.lognot x) in
+        Some { bitstring ; phase }
+      end
+
+let single_excitation_reference h p spindet =
+      creation p @@ annihilation h @@ spindet
+
+let single_excitation h p =
+  single_excitation_reference h p
+
+
+let double_excitation_reference h' p' h p spindet =
+      creation p' @@ creation p @@ annihilation h @@ annihilation h' @@ spindet
+
+let double_excitation h' p' h p =
+  double_excitation_reference h' p' h p
+
+
+let excitation_level t t' =
+  Bitstring.hamdist (bitstring t) (bitstring t') / 2
+
+let holes_of t t' =
+  Bitstring.logand (bitstring t) (Bitstring.logxor (bitstring t) (bitstring t'))
+  |> Bitstring.to_list
+
+let particles_of t t' =
+  Bitstring.logand (bitstring t') (Bitstring.logxor (bitstring t) (bitstring t'))
+  |> Bitstring.to_list
+
+let holes_particles_of t t' =
+  let x = Bitstring.logxor (bitstring t) (bitstring t') in
+  let holes = Bitstring.logand (bitstring t) x |> Bitstring.to_list
+  and particles = Bitstring.logand (bitstring t') x |> Bitstring.to_list
+  in
+  List.rev_map2 (fun h p -> (h,p)) holes particles
+  |> List.rev
+
+
+let set_phase p = function
+  | Some t -> Some { t with phase = p }
+  | None -> None
+
+let negate_phase = function
+  | Some t -> Some { t with phase = Phase.neg t.phase }
+  | None -> None
+
+
+let of_bitstring ?(phase=Phase.Pos) bitstring =
+  Some { bitstring ; phase }
+
+let of_list n l =
+  List.rev l
+  |> List.fold_left (fun accu p -> creation p accu) (vac n)
+
+
+let to_list = function
+  | None -> []
+  | Some spindet ->
+    let rec aux accu z =
+      if not (Bitstring.is_zero z) then
+        let element = ((Bitstring.trailing_zeros z)+1) in
+        (aux [@tailcall]) (element::accu) (Bitstring.logand z (Bitstring.minus_one z) )
+      else List.rev accu
+    in aux [] spindet.bitstring
+
+let to_array t =
+  to_list t
+  |> Array.of_list
+
+let n_electrons = function
+  | Some t -> Bitstring.popcount t.bitstring
+  | None -> 0
+
+
+let pp ppf = function
+  | None -> Format.fprintf ppf "@[<h>None@]"
+  | Some s ->
+      Format.fprintf ppf "@[<h>%a %a@]" Phase.pp s.phase Bitstring.pp
+      s.bitstring
+
+

ci/lib/spindeterminant.mli

@@ -0,0 +1,91 @@
+(**
+A spin-determinant is one of the two determinants in the Waller-Hartree
+double determinant representation of a Slater determinant. It is represented
+as a bit string and a phase factor.
+*)
+
+open Common
+
+type t
+type hole = int
+type particle = int
+
+(** {1 Accessors}. *)
+
+val phase : t -> Phase.t
+(** Phase factor.
+    @raise Invalid_argument if the spin-determinant is [None].
+*)
+
+val set_phase : Phase.t -> t -> t
+(** Returns a spin-determinant with the phase set to [p]. *)
+
+
+val bitstring : t -> Bitstring.t
+(** Bit string.
+    @raise Invalid_argument if the spin-determinant is [None].
+*)
+
+val is_none : t -> bool
+(** Tests if a spin-determinant is [None]. *)
+
+val negate_phase : t -> t
+(** Returns a spin-determinant with the phase reversed. *)
+
+
+(** {1 Second quantization operators} *)
+
+val vac : int -> t
+(** Vacuum state, [vac = Some ]{% $|\rangle$ %}. The integer parameter contains the
+   number of orbitals in the basis set. *)
+
+val creation : particle -> t -> t
+(** [creation p] is the creation operator {% $a^\dagger_p$ %}. *)
+
+val annihilation : hole -> t -> t
+(** [annihilation h] is the annihilation operator {% $a_h$ %}. *)
+
+val single_excitation : hole -> particle -> t -> t
+(** Single excitation operator {% $T_h^p = a^\dagger_p a_h$ %}. *)
+
+val double_excitation : hole -> particle -> hole -> particle -> t -> t
+(** Double excitation operator {% $T_{hh'}^{pp'} = a^\dagger_p a^\dagger_{p'} a_{h'} a_h$ %}. *)
+
+val excitation_level : t -> t -> int
+(** Returns the excitation level between two spin-determinants. *)
+
+val holes_of : t -> t -> int list
+(** Returns the list of holes in the excitation from one determinant to another. *)
+
+val particles_of : t -> t -> int list
+(** Returns the list of particles in the excitation from one determinant to another. *)
+
+val holes_particles_of : t -> t -> (int*int) list
+(** Returns the list of pairs of holes/particles in the excitation from one determinant to
+another. *)
+
+val n_electrons : t -> int
+(** Returns the number of electrons in the determinant. *)
+
+
+(** {1 Creation} *)
+
+val of_bitstring : ?phase:Phase.t -> Bitstring.t -> t
+(** Creates from a bitstring and an optional phase.*)
+
+val of_list : int -> 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]
+    spin-determinant is returned.
+    The first integer is the size of the MO basis set. *)
+
+val to_list : t -> int list
+(** Transforms a spin-determinant into a list of orbital indices. *)
+
+val to_array : t -> int array
+(** Transforms a spin-determinant into an array of orbital indices. *)
+
+(** {1 Printers}. *)
+
+val pp : Format.formatter -> t -> unit
+

ci/lib/spindeterminant_space.ml

@@ -0,0 +1,75 @@
+open Common
+
+type t =
+  {
+    elec_num          : int;
+    mo_basis          : Mo.Basis.t;
+    mo_class          : Mo.Class.t;
+    spin_determinants : Spindeterminant.t array;
+  }
+
+
+let spin_determinants t = t.spin_determinants
+let elec_num t = t.elec_num
+let mo_basis t = t.mo_basis
+let mo_class t = t.mo_class
+let size     t = Array.length t.spin_determinants
+
+let fci_of_mo_basis ~frozen_core mo_basis elec_num =
+  let mo_num = Mo.Basis.size mo_basis in
+  let mo_class = Mo.Class.fci ~frozen_core mo_basis in
+  let m l = 
+    List.fold_left (fun accu i -> let j = i-1 in
+      Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
+                   ) (Bitstring.zero mo_num) l
+  in
+  let occ_mask    = m (Mo.Class.core_mos    mo_class)
+  and active_mask = m (Mo.Class.active_mos  mo_class)
+  in
+  let neg_active_mask = Bitstring.lognot active_mask in
+  (* Here we generate the FCI space and filter out unwanted determinants
+     with excitations involving the core electrons. This should be improved. *)
+  let spin_determinants = 
+    Bitstring.permutations elec_num mo_num 
+    |> List.filter (fun b -> Bitstring.logand neg_active_mask b = occ_mask)
+    |> Array.of_list
+    |> Array.map (fun b -> Spindeterminant.of_bitstring b)
+  in
+  { elec_num ; mo_basis ; mo_class ; spin_determinants }
+
+
+let cas_of_mo_basis mo_basis ~frozen_core elec_num n m =
+  let mo_num = Mo.Basis.size mo_basis in
+  let mo_class = Mo.Class.cas_sd ~frozen_core mo_basis n m in
+  let m l = 
+    List.fold_left (fun accu i -> let j = i-1 in
+      Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
+                   ) (Bitstring.zero mo_num) l
+  in
+  let active_mask   = m (Mo.Class.active_mos   mo_class) in
+  let occ_mask      = m (Mo.Class.core_mos     mo_class) in
+  let inactive_mask = m (Mo.Class.inactive_mos mo_class) in
+  let occ_mask = Bitstring.logor occ_mask inactive_mask in
+  let neg_active_mask = Bitstring.lognot active_mask in
+  (* Here we generate the FCI space and filter out all the unwanted determinants.
+     This should be improved. *)
+  let spin_determinants = 
+    Bitstring.permutations elec_num mo_num 
+    |> List.filter (fun b -> Bitstring.logand neg_active_mask b = occ_mask)
+    |> Array.of_list
+    |> Array.map (fun b -> Spindeterminant.of_bitstring b)
+  in
+  { elec_num ; mo_basis ; mo_class ; spin_determinants }
+
+
+
+let pp ppf t =
+  Format.fprintf ppf "@[<v 2> [";
+  Array.iteri (fun i d ->
+    Format.fprintf ppf "@[<v>@[%8d@] @[%a@]@]@;" i Spindeterminant.pp d)
+    (spin_determinants t) ;
+  Format.fprintf ppf "]@]"
+
+
+
+

ci/lib/spindeterminant_space.mli

@@ -0,0 +1,46 @@
+(**
+The space built with determinants made with same-spin spinorbitals.
+*)
+
+type t
+
+(** {1 Accessors} *)
+
+val size : t -> int
+(** Number of determinants in the space. *)
+
+val spin_determinants : t -> Spindeterminant.t array
+(** All the spin-determinants belonging to the space. *)
+
+val elec_num : t -> int
+(** Number of (same-spin) electrons occupying the 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. *)
+
+
+(** {1 Creation} *)
+
+val fci_of_mo_basis : frozen_core:Mo.Frozen_core.t -> Mo.Basis.t -> int -> t
+(** Create a space of all possible ways to put [n_elec-ncore] electrons in the
+    [Active] MOs.  All other MOs are untouched.
+*)
+
+val cas_of_mo_basis : Mo.Basis.t -> frozen_core:Mo.Frozen_core.t-> int -> int -> int -> t
+(** [cas_of_mo_basis mo_basis n_elec n m] creates a CAS(n,m) space of 
+    [Active] MOs. The unoccupied MOs are [Virtual], and the occupied MOs
+    are [Core] and [Inactive].
+*)
+
+
+(** {2 Printing} *)
+
+val pp : Format.formatter -> t -> unit
+
+
+
+
+

ci/test/determinant.ml

@@ -0,0 +1,111 @@
+open Common
+open Ci
+
+let tests =
+
+  let test_creation () =
+      let l_a = [ 1 ; 2 ; 3 ; 5 ; 64 ]
+      and l_b = [ 2 ; 3 ; 5 ; 65 ] in
+      let det = Determinant.of_lists 66 l_a l_b in
+      let z_a = Determinant.alfa det
+      and z_b = Determinant.beta  det in
+      Alcotest.(check (list int )) "alfa" (Spindeterminant.to_list z_a) l_a;
+      Alcotest.(check (list int )) "beta"  (Spindeterminant.to_list z_b) l_b;
+      Alcotest.(check bool) "phase" (Determinant.phase det = Phase.Pos)  true;
+  in
+
+  let test_phase () =
+      let l_a = [ 1 ; 2 ; 3 ; 64 ; 5 ]
+      and l_b = [ 2 ; 3 ; 5 ; 65 ] in
+      let det = Determinant.of_lists 66 l_a l_b in
+      Alcotest.(check bool) "phase" (Determinant.phase det = Phase.Neg)  true;
+      let l_a = [ 1 ; 2 ; 3 ; 64 ; 5 ]
+      and l_b = [ 3 ; 2 ; 5 ; 65 ] in
+      let det = Determinant.of_lists 66 l_a l_b in
+      Alcotest.(check bool) "phase" (Determinant.phase det = Phase.Pos)  true;
+      let l_a = [ 1 ; 3 ; 2 ; 64 ; 5 ]
+      and l_b = [ 3 ; 2 ; 5 ; 65 ] in
+      let det = Determinant.of_lists 66 l_a l_b in
+      Alcotest.(check bool) "phase" (Determinant.phase det = Phase.Neg)  true;
+      let l_a = [ 1 ; 3 ; 2 ; 64 ; 5 ]
+      and l_b = [ 3 ; 2 ; 65 ; 5 ] in
+      let det = Determinant.of_lists 66 l_a l_b in
+      Alcotest.(check bool) "phase" (Determinant.phase det = Phase.Pos)  true;
+  in
+
+  let test_operators () =
+      let det =
+        let open Determinant in
+        let open Spin in
+        creation Alfa 1 @@ creation Alfa 3 @@ creation Alfa 2 @@ creation Alfa 5 @@ 
+        creation Beta 1 @@ creation Beta 3 @@ creation Beta 4 @@ creation Beta 5 @@ vac 10
+      in
+      Alcotest.(check bool) "creation 1" true
+        (det = Determinant.of_lists 10 [ 1 ; 3 ; 2 ; 5 ] [1 ; 3 ; 4 ; 5 ] );
+
+      let det' =
+        Determinant.single_excitation Spin.Alfa 3 6 det
+      in
+      Alcotest.(check bool) "single_exc 1" true
+        (det' = Determinant.of_lists 10 [ 1 ; 6 ; 2 ; 5 ] [1 ; 3 ; 4 ; 5 ] );
+
+      let det' =
+        Determinant.single_excitation Spin.Beta  3 6 det
+      in
+      Alcotest.(check bool) "single_exc 2" true
+        (det' = Determinant.of_lists 10 [ 1 ; 3 ; 2 ; 5 ] [1 ; 6 ; 4 ; 5 ] );
+
+      let det' =
+        Determinant.single_excitation Spin.Alfa 4 6 det
+      in
+      Alcotest.(check bool) "single_exc 3" true (Determinant.is_none det');
+
+      let det' =
+        Determinant.single_excitation Spin.Beta  1 5 det
+      in
+      Alcotest.(check bool) "single_exc 4" true (Determinant.is_none det');
+
+      let det' =
+        Determinant.double_excitation Spin.Alfa 3 6 Spin.Alfa 2 7 det
+      in
+      let det'' = Determinant.of_lists 10 [ 1 ; 6 ; 7 ; 5 ] [1 ; 3 ; 4 ; 5 ] in
+      Alcotest.(check bool) "double_exc 1" true (det' = det'');
+
+      let det' =
+        Determinant.double_excitation Spin.Beta  3 6 Spin.Beta 5 7  det
+      in
+      Alcotest.(check bool) "double_exc 2" true
+        (det' = Determinant.of_lists 10 [ 1 ; 3 ; 2 ; 5 ] [1 ; 6 ; 4 ; 7 ] );
+
+      let det' =
+        Determinant.double_excitation Spin.Alfa 3 6 Spin.Beta  5 7 det
+      in
+      Alcotest.(check bool) "double_exc 3" true
+        (det' = Determinant.of_lists 10 [ 1 ; 6 ; 2 ; 5 ] [1 ; 3 ; 4 ; 7 ] );
+
+      let det' =
+        Determinant.double_excitation Spin.Beta 5 7 Spin.Alfa 3 6 det
+      in
+      Alcotest.(check bool) "double_exc 4" true
+        (det' = Determinant.of_lists 10 [ 1 ; 6 ; 2 ; 5 ] [1 ; 3 ; 4 ; 7 ] );
+
+      let det' =
+        Determinant.double_excitation Spin.Alfa 4 6 Spin.Alfa 2 7 det
+      in
+      Alcotest.(check bool) "double_exc 5" true (Determinant.is_none det');
+
+      let det' =
+        Determinant.double_excitation Spin.Beta  1 5 Spin.Alfa 2 7  det
+      in
+      Alcotest.(check bool) "double_exc 6" true (Determinant.is_none det');
+
+  in
+
+  [
+    "Creation", `Quick, test_creation;
+    "Phase",    `Quick, test_phase;
+    "Operators",`Quick, test_operators;
+  ]
+
+
+  

ci/test/dune

@@ -0,0 +1,9 @@
+(library
+ (name test_ci)
+ (synopsis "Test for ci library")
+ (libraries
+  alcotest
+  trexio
+  qcaml.ci
+ )
+)

ci/test/excitation.ml

@@ -0,0 +1,52 @@
+open Common
+open Ci
+
+let tests =
+
+  (*
+  let test_id () =
+    let l_a = [ 1 ; 2 ; 3 ; 5 ; 64 ]
+    and l_b = [ 2 ; 3 ; 5 ; 65 ] in
+    let det1 = Determinant.of_lists l_a l_b in
+    let det2 = Determinant.negate_phase det1 in
+  in
+  *)
+
+  let test_single () =
+    let l_a = [ 1 ; 2 ; 3 ; 5 ; 64 ]
+    and l_b = [ 2 ; 3 ; 5 ; 65 ] in
+    let det1 = Determinant.of_lists 66 l_a l_b in
+    let det2 = Determinant.single_excitation Spin.Alfa 3 7 det1 in
+    let t = Excitation.single_of_det det1 det2 in
+    Alcotest.(check bool) "single 1" true (t = Single (Phase.Pos, { hole=3 ; particle=7 ; spin=Spin.Alfa}) );
+    let det2 =
+      Determinant.single_excitation Spin.Alfa 2 7 det1
+      |> Determinant.negate_phase
+    in
+    let t = Excitation.single_of_det det1 det2 in
+    Alcotest.(check bool) "single 2" true (t = Single (Phase.Neg, { hole=2 ; particle=7 ; spin=Spin.Alfa }) );
+    let det2 = Determinant.single_excitation Spin.Beta 2 7 det1 in
+    let t = Excitation.single_of_det det1 det2 in
+    Alcotest.(check bool) "single 3" true (t = Single (Phase.Pos, { hole=2 ; particle=7 ; spin=Spin.Beta}) );
+    let det2 = Determinant.single_excitation Spin.Beta 3 256 det1 in
+    let t = Excitation.single_of_det det1 det2 in
+    Alcotest.(check bool) "single 4" true (t = Single (Phase.Pos, { hole=3 ; particle=256 ; spin=Spin.Beta}) );
+  in
+
+  let test_double () =
+    let l_a = [ 1 ; 2 ; 3 ; 5 ; 64 ]
+    and l_b = [ 2 ; 3 ; 5 ; 65 ] in
+    let det1 = Determinant.of_lists 66 l_a l_b in
+    let det2 = Determinant.double_excitation Spin.Alfa 3 7 Spin.Alfa 2 6 det1 in
+    let t = Excitation.double_of_det det1 det2 in
+    Alcotest.(check bool) "double 1" true
+      (t = Double (Phase.Neg,
+        { hole=2 ; particle=7 ; spin=Spin.Alfa},
+        { hole=3 ; particle=6 ; spin=Spin.Alfa}));
+  in
+  [
+    "Single", `Quick, test_single;
+    "Double", `Quick, test_double;
+  ]
+
+

ci/test/spindeterminant.ml

@@ -0,0 +1,102 @@
+open Alcotest
+open Ci
+open Ci.Spindeterminant
+
+let tests =
+
+  let test_creation () =
+    let l_a = [ 1 ; 2 ; 3 ; 5 ] in
+    let det = of_list 10 l_a in
+    check (list int ) "bitstring 1" l_a (to_list det);
+    check bool "phase 2" true (phase det = Phase.Pos);
+    let l_b = [ 1 ; 3 ; 2 ; 5 ] in
+    let det = of_list 10 l_b in
+    check (list int ) "bitstring 2" l_a (to_list det);
+    check bool "phase 2" true (phase det = Phase.Neg);
+  in
+
+  let test_a_operators () =
+    let det =
+      creation 5 @@ creation 2 @@ creation 2 @@ creation 1 @@ (vac 10)
+    in
+    check bool "none 1" true (is_none det);
+
+    let det =
+      creation 5 @@ creation 3 @@ creation 2 @@ creation 1 @@ (vac 10)
+    in
+    let l_a = [ 1 ; 2 ; 3 ; 5 ] in
+    check (list int) "bitstring 1" l_a (to_list det);
+    check bool "phase 1" true (phase det = Phase.Pos);
+
+    let det =
+      creation 1 @@ creation 3 @@ creation 2 @@ creation 5 @@ (vac 10)
+    in
+    check (list int) "bitstring 2" l_a (to_list det);
+    check bool "phase 2" true (phase det = Phase.Neg);
+
+    let l_b = [ 1 ; 3 ; 2 ; 5 ] in
+    let det = of_list 10 l_b in
+    check (list int ) "bitstring 3" l_a (to_list det);
+    check bool "phase 3" true (phase det = Phase.Neg);
+
+    check bool "none 1" true (annihilation 4 det |> is_none);
+
+    let det =
+      annihilation 1 det
+    in
+    check (list int ) "bitstring 4" (List.tl l_a) (to_list det);
+    check bool "phase 4" true (phase det = Phase.Neg);
+
+    let det =
+      annihilation 3 det
+    in
+    check (list int ) "bitstring 5" [ 2 ; 5 ] (to_list det);
+    check bool "phase 5" true (phase det = Phase.Pos);
+
+    let det =
+      annihilation 5 @@ annihilation 2 det
+    in
+    check (list int ) "bitstring 6" [] (to_list det);
+    check bool "phase 6" true (phase det = Phase.Pos);
+
+  in
+
+  let test_exc_operators () =
+    let l_a = [ 1 ; 2 ; 3 ; 5 ] in
+    let det = of_list 10 l_a in
+    let l_b = [ 1 ; 7 ; 3 ; 5 ] in
+    let det2 = of_list 10 l_b in
+    Format.printf "%a@." pp det;
+    Format.printf "%a@." pp det2;
+    Format.printf "%a@." pp (Spindeterminant.single_excitation 2 7 det);
+    check bool "single 1" true (Spindeterminant.single_excitation 2 7 det = det2);
+    check bool "single 2" true (Spindeterminant.single_excitation 4 7 det |> is_none);
+
+    let l_c = [ 1 ; 7 ; 6 ; 5 ] in
+    let det3 = of_list 10 l_c in
+    check bool "double 1" true (double_excitation 2 7 3 6 det = det3);
+    check bool "double 2" true (double_excitation 4 7 3 6 det |> is_none);
+  in
+
+  let test_exc_spindet () =
+    let l_a = [ 1 ; 2 ; 3 ; 5 ] in
+    let det = of_list 10 l_a in
+    let l_b = [ 1 ; 7 ; 3 ; 5 ] in
+    let det2 = of_list 10 l_b in
+    check int "single" 1 (excitation_level det det2);
+    check (list int) "holes" [2] (holes_of det det2);
+    check (list int) "particles" [7] (particles_of det det2);
+    let l_b = [ 1 ; 7 ; 3 ; 6 ] in
+    let det2 = of_list 10 l_b in
+    check int "double" 2 (excitation_level det det2);
+    check (list int) "holes" [2 ; 5] (holes_of det det2);
+    check (list int) "particles" [6 ; 7] (particles_of det det2);
+  in
+  [
+    "Creation", `Quick, test_creation;
+    "Creation/Annihilation Operators", `Quick, test_a_operators;
+    "Excitation Operators", `Quick, test_exc_operators;
+    "Excitation of spindet", `Quick, test_exc_spindet;
+  ]
+
+

examples/ex_trexio_integrals.org

@@ -0,0 +1,204 @@
+#+TITLE: Integrals
+
+#+PROPERTY
+
+In this example, we write a program that reads the geometry of a
+molecule in =xyz= format and a Gaussian atomic basis set in GAMESS
+format. The output is written in a trexio file.
+
+* Header 
+
+#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+module Command_line = Qcaml.Common.Command_line
+module Util = Qcaml.Common.Util
+open Qcaml.Linear_algebra
+#+END_SRC
+
+#+RESULTS:
+: module Command_line = Qcaml.Common.Command_line
+: module Util = Qcaml.Common.Util
+
+#+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let () =
+#+END_SRC
+
+* Command-line arguments
+
+  We use the =Command_line= module to define the following possible
+  arguments:
+  - =-b --basis= : The name of the file containing the basis set
+  - =-x --xyz=   : The name of the file containing the atomic coordinates
+  - =-u --range-separation= : The value of $\mu$, the range-separation
+    parameter in range-separated DFT. If this option is not present,
+    no output file will be generated for the range-separated integrals.
+
+** Definition
+
+   #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let open Command_line in
+begin
+  set_header_doc (Sys.argv.(0));
+  set_description_doc "Computes the one- and two-electron integrals on the Gaussian atomic basis set.";
+  set_specs
+    [ { short='b' ; long="basis" ; opt=Mandatory;
+	arg=With_arg "<string>";
+	doc="Name of the file containing the basis set"; } ;
+
+      { short='x' ; long="xyz" ; opt=Mandatory;
+	arg=With_arg "<string>";
+	doc="Name of the file containing the nuclear coordinates in xyz format"; } ;
+
+      { short='u' ; long="range-separation" ; opt=Optional;
+	arg=With_arg "<float>";
+	doc="Range-separation parameter."; } ;
+    ]
+end;
+   #+END_SRC
+
+** Interpretation
+
+   #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let basis_file  = Util.of_some @@ Command_line.get "basis" in
+let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
+
+let range_separation =
+  match Command_line.get "range-separation" with
+  | None -> None
+  | Some mu -> Some (float_of_string mu)
+in
+
+let operators =
+  match range_separation with
+  | None -> []
+  | Some mu -> [ Qcaml.Operators.Operator.of_range_separation mu ]
+in
+   #+END_SRC
+
+* Computation
+  We first read the =xyz= file to create a molecule:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let nuclei =
+  Qcaml.Particles.Nuclei.of_xyz_file nuclei_file
+in
+  #+END_SRC
+
+  Then we create an Gaussian AO basis using the atomic coordinates,
+  and we optionally introduce the range-separation parameter via the
+ =operators=:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let ao_basis =
+  Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian
+    ~operators ~cartesian:true ~nuclei basis_file
+in
+  #+END_SRC
+
+  We compute the required one-electron integrals:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let overlap   = Qcaml.Ao.Basis.overlap   ao_basis in
+let eN_ints   = Qcaml.Ao.Basis.eN_ints   ao_basis in
+let kin_ints  = Qcaml.Ao.Basis.kin_ints  ao_basis in
+let multipole = Qcaml.Ao.Basis.multipole ao_basis in
+let x_mat     = multipole "x" in
+let y_mat     = multipole "y" in
+let z_mat     = multipole "z" in
+  #+END_SRC
+
+  and the two-electron integrals (1/r and long range):
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+let ee_ints   = Qcaml.Ao.Basis.ee_ints   ao_basis in
+let lr_ints   = 
+  match range_separation with
+  | Some _mu -> Some (Qcaml.Ao.Basis.ee_lr_ints ao_basis)
+  | None -> None
+in
+  #+END_SRC
+
+* Output
+
+  We write the one-electron integrals:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+Matrix.to_file ~filename:"overlap.dat" ~sym:true overlap;
+Matrix.to_file ~filename:"eN.dat" ~sym:true eN_ints;
+Matrix.to_file ~filename:"kinetic.dat" ~sym:true kin_ints;
+Matrix.to_file ~filename:"x.dat" ~sym:true x_mat;
+Matrix.to_file ~filename:"y.dat" ~sym:true y_mat;
+Matrix.to_file ~filename:"z.dat" ~sym:true z_mat;
+  #+END_SRC
+
+  and the the two-electron integrals:
+  #+BEGIN_SRC ocaml :comments link :exports code :tangle ex_integrals.ml  
+Four_idx_storage.to_file ~filename:"eri.dat" ee_ints;
+match lr_ints with
+| Some integrals -> Four_idx_storage.to_file ~filename:"eri_lr.dat" integrals;
+| None -> ()
+  #+END_SRC
+  
+
+* Interactive test                                                 :noexport:
+  #+begin_src ocaml :results drawer :session :cache no  :exports none
+#require "qcaml.top" ;;
+#require "trexio" ;;
+open Qcaml ;;
+  #+end_src
+
+  #+RESULTS:
+  :results:
+  :end:
+
+  #+begin_src ocaml :results drawer :session :cache no :exports none
+let nuclei_file = "/dev/shm/f2.xyz"
+let nuclei =
+  Qcaml.Particles.Nuclei.of_xyz_file nuclei_file
+  #+end_src
+
+  #+RESULTS:
+  :results:
+  val nuclei_file : string = "/dev/shm/f2.xyz"
+  val nuclei : Qcaml.Particles.Nuclei.t =
+
+                  Nuclear Coordinates (Angstrom)
+                  ------------------------------
+
+  -----------------------------------------------------------------------
+     Center   Atomic   Element          Coordinates (Angstroms)
+              Number                 X             Y             Z
+  -----------------------------------------------------------------------
+       1        9        F       0.000000      0.000000      0.000000
+       2        9        F       0.000000      0.000000      1.411900
+  -----------------------------------------------------------------------
+  :end:
+
+  #+begin_src ocaml :results drawer :session :cache no :exports none
+let trexio_file = Trexio.open_file "test.trexio" 'w' Trexio.HDF5
+  #+end_src
+
+  #+RESULTS:
+  :results:
+  val trexio_file : Trexio.trexio_file = <abstr>
+  :end:
+
+  #+begin_src ocaml :results drawer :session :cache no :exports none
+Qcaml.Particles.Nuclei.to_trexio trexio_file nuclei;;
+Trexio.close_file trexio_file;;
+  #+end_src
+
+  #+RESULTS:
+  :results:
+  - : unit = ()
+  :end:
+
+  #+begin_src ocaml :results drawer :session :cache no :exports none
+let basis_file = "/home/scemama/qp2/data/basis/cc-pvdz";; 
+let ao_basis =
+  Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file
+  #+end_src
+
+  #+RESULTS:
+  :results:
+  val ao_basis : Qcaml.Ao.Basis.t = Gaussian Basis, spherical, 30 AOs
+  :end:
+
+
+  #+begin_src ocaml :results drawer :session :cache no :exports none
+Trexio.close_file trexio_file;;
+  #+end_src

qcaml/lib/dune

@@ -7,6 +7,7 @@
   )
  (libraries
  qcaml.ao
+ qcaml.ci
  qcaml.common
  qcaml.gaussian
  qcaml.gaussian_integrals

qcaml/lib/qcaml.ml

@@ -1,5 +1,5 @@
-(* Auto-generated by qcaml/README.org *)
 module Ao = Ao
+module Ci = Ci
 module Common = Common
 module Gaussian = Gaussian
 module Gaussian_integrals = Gaussian_integrals

test/dune

@@ -9,6 +9,7 @@
   test_gaussian_integrals
   test_ao
   test_mo
+  test_ci
   test_perturbation
 ))
 

test/run_tests.ml

@@ -14,6 +14,9 @@ let test_suites: unit Alcotest.test list = [
   "Mo.Guess", Test_mo.Guess.tests;
   "Mo.Hartree_Fock", Test_mo.Hf.tests;
   "Perturbation.Mp2", Test_perturbation.Mp2.tests;
+  "Ci.Spindeterminant", Test_ci.Spindeterminant.tests;
+  "Ci.Determinant", Test_ci.Determinant.tests;
+  "Ci.Excitation", Test_ci.Excitation.tests;
 ]
 
 let () = Alcotest.run "QCaml" test_suites