Determinant space OK

Basis/PrimitiveShellPair.ml

@@ -3,15 +3,15 @@ open Constants
 
 
 type t = {
-  exponent        : float;              (* alpha + beta *)
-  exponent_inv    : float;              (* 1/(alpha + beta) *)
-  center          : Coordinate.t;       (* P = (alpha * A + beta B)/(alpha+beta) *)
-  center_minus_a  : Coordinate.t;       (* P - A *)
-  a_minus_b       : Coordinate.t;       (* A - B *)
-  a_minus_b_sq    : float;              (* |A-B|^2 *)
+  exponent        : float;              (* {% $\alpha + \beta$ %} *)
+  exponent_inv    : float;              (* {% $1/(\alpha + \beta)$ %} *)
+  center          : Coordinate.t;       (* {% $P = (\alpha A + \beta B)/(\alpha+\beta)$ %} *)
+  center_minus_a  : Coordinate.t;       (* {% $P - A$ %} *)
+  a_minus_b       : Coordinate.t;       (* {% $A - B$ %} *)
+  a_minus_b_sq    : float;              (* {% $|A-B|^2$ %} *)
   norm_scales     : float array lazy_t;
-  normalization   : float;              (* norm_coef_a * norm_coef_b * g, with
-                     g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *)
+  normalization   : float;              (* [norm_coef_a * norm_coef_b * g], with
+                     {% $g = (\pi/(\alpha+\beta))^(3/2) \exp (-|A-B|^2 \alpha\beta/(\alpha+\beta))$ %} *)
   ang_mom         : AngularMomentum.t; 
   shell_a         : PrimitiveShell.t;
   shell_b         : PrimitiveShell.t;
@@ -69,7 +69,7 @@ let create_make_of p_a p_b =
       Ps.normalization p_a
     in
 
-    let alpha_a = 
+    let alfa_a = 
       Co.( Ps.exponent p_a |. Ps.center p_a )
     in
 
@@ -101,7 +101,7 @@ let create_make_of p_a p_b =
           in
 
           let center =
-            Co.(exponent_inv |. (alpha_a |+ beta_b))
+            Co.(exponent_inv |. (alfa_a |+ beta_b))
           in
 
           let center_minus_a = 

CI/Determinant.ml

@@ -76,11 +76,11 @@ let double_excitation spin h p spin' h' p' t =
 
 
 
-let pp_det ppf t =
-  Format.fprintf ppf "@[<v>@[phase:%a@]@;@[a:%a@];@[b:%a@]@]@."
+let pp_det n ppf t =
+  Format.fprintf ppf "@[<v>@[phase:%a@]@;@[a:%a@]@;@[b:%a@]@]@."
     Phase.pp_phase (phase t)
-    Spindeterminant.pp_spindet t.alfa
-    Spindeterminant.pp_spindet t.beta
+    (Spindeterminant.pp_spindet n) t.alfa
+    (Spindeterminant.pp_spindet n) t.beta
 
 
 

CI/Determinant.mli

@@ -62,7 +62,8 @@ val negate_phase : t -> t
 
 (** {1 Printers} *)
 
-val pp_det : Format.formatter -> t -> unit
+val pp_det : int -> Format.formatter -> t -> unit
+(** First [int] is the number of MOs to print. *)
 
 (** {1 Unit tests} *)
 

CI/Determinant_space.ml

@@ -0,0 +1,48 @@
+type t = 
+{
+  n_alfa       : int ;
+  n_beta       : int ;
+  mo_class     : MOClass.t ;
+  mo_basis     : MOBasis.t ;
+  determinants : Determinant.t array;
+}
+
+module Ss = Spindeterminant_space
+
+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 determinants t = t.determinants
+
+let fci_of_mo_basis ?(frozen_core=true) mo_basis =
+  let s = MOBasis.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 =
+    Ss.fci_of_mo_basis ~frozen_core mo_basis n_alfa
+  and det_b =
+    Ss.fci_of_mo_basis ~frozen_core mo_basis n_beta
+  in
+  let mo_class = Ss.mo_class det_a in
+  let determinants =
+    let a = Ss.spin_determinants det_a
+    and b = Ss.spin_determinants det_b
+    in
+    let sb = Ss.size det_b in
+    Array.to_list a
+    |> List.map (fun a_i ->
+        Array.init sb (fun j ->
+            Determinant.of_spindeterminants a_i b.(j)) )
+    |> Array.concat
+  in
+  { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants }
+
+
+let pp_det_space ppf t =
+  Format.fprintf ppf "@[<v 2>[ ";
+  Array.iteri (fun i d -> Format.fprintf ppf "@[<v>@[%8d@]@;@[%a@]@]@;" i
+                    (Determinant.pp_det (MOBasis.size (mo_basis t))) d) t.determinants ;
+  Format.fprintf ppf "]@]"
+

CI/Determinant_space.mli

@@ -0,0 +1,32 @@
+(**
+The determinant space in which we solve the Schrodinger equation.
+*)
+
+type t
+
+(** {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 -> MOClass.t
+(** The MO classes used to generate the space. *)
+
+val mo_basis : t -> MOBasis.t
+(** The MO basis on which the determinants are expanded. *)
+
+val determinants : t -> Determinant.t array
+(** All the determinants belonging to the space. *)
+
+val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.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.
+*)
+
+(** {2 Printing} *)
+
+val pp_det_space : Format.formatter -> t -> unit

CI/Spindeterminant.ml

@@ -116,10 +116,10 @@ let rec to_list = function
     in aux [] spindet.bitstring
 
 
-let pp_spindet ppf = function
+let pp_spindet n ppf = function
   | None -> Format.fprintf ppf "@[<h>None@]"
   | Some s -> 
-      Format.fprintf ppf "@[<h>%a %a@]" Phase.pp_phase s.phase Util.pp_bitstring s.bitstring
+      Format.fprintf ppf "@[<h>%a %a@]" Phase.pp_phase s.phase (Util.pp_bitstring n) s.bitstring
 
 
 

CI/Spindeterminant.mli

@@ -69,7 +69,8 @@ val to_list : t -> int list
 
 (** {1 Printers}. *)
 
-val pp_spindet : Format.formatter -> t -> unit
+val pp_spindet : int -> Format.formatter -> t -> unit
+(** First [int] is the number of MOs to print *)
 
 
 (** {1 Unit testing} *)

CI/Spindeterminant_space.ml

@@ -1,48 +1,42 @@
 type t =
-{
-  elec_num : int;
-  mo_basis : MOBasis.t;
-  mo_class : MOClass.t;
-  spindets : Spindeterminant.t array;
-}
+  {
+    elec_num          : int;
+    mo_basis          : MOBasis.t;
+    mo_class          : MOClass.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=true) mo_basis elec_num =
   let mo_num = MOBasis.size mo_basis in
-  let ncore  = (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 in
-  let mo_class =
-    MOClass.of_list @@
-    if frozen_core then
-      List.concat [
-        Util.list_range 1 ncore
-        |> List.map (fun i -> MOClass.Core i) ;
-        Util.list_range (ncore+1) mo_num
-        |> List.map (fun i -> MOClass.Active i) 
-      ]
-    else
-      Util.list_range 1 mo_num
-      |> List.map (fun i -> MOClass.Active i)
-  in
+  let mo_class = MOClass.fci ~frozen_core mo_basis in
   let m l = 
     List.fold_left (fun accu i -> let j = i-1 in Z.(logor accu (shift_left one j)) 
-      ) Z.zero l
+                   ) Z.zero l
   in
   let occ_mask    = m (MOClass.core_mos    mo_class)
   and active_mask = m (MOClass.active_mos  mo_class)
   in
   let neg_active_mask = Z.lognot active_mask in
-  let spindets = 
+  let spin_determinants = 
     Util.bit_permtutations elec_num mo_num 
     |> List.filter (fun b -> Z.logand neg_active_mask b = occ_mask)
     |> List.map (fun b -> Spindeterminant.of_bitstring b)
     |> Array.of_list
   in
-  { elec_num ; mo_basis ; mo_class ; spindets }
+  { elec_num ; mo_basis ; mo_class ; spin_determinants }
 
 
 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 "@[<v 2>[ ";
+  Array.iteri (fun i d -> Format.fprintf ppf "@[<v>@[%8d@] @[%a@]@]@;" i
+                  (Spindeterminant.pp_spindet (MOBasis.size (mo_basis t))) d) (spin_determinants t) ;
   Format.fprintf ppf "]@]"
 
 

CI/Spindeterminant_space.mli

@@ -1,11 +1,28 @@
-type t =
-{
-  elec_num : int                     ; (** Number of electrons *)
-  mo_basis : MOBasis.t               ; (** MO basis on which the space is built *)
-  mo_class : MOClass.t               ; (** CI Classes of the MOs *)
-  spindets : Spindeterminant.t array ; (** Spin-determinants belonging to the space *)
-}
+(**
+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 -> MOClass.t
+(** The MO classes used to generate the space. *)
+
+val mo_basis : t -> MOBasis.t
+(** The MO basis on which the determinants are expanded. *)
+
+
+(** {1 Creation} *)
 
 val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t
 (** Create a space of all possible ways to put [n_elec] electrons in the [Active] MOs.

MOBasis/MOClass.ml

@@ -55,3 +55,19 @@ let deleted_mos t =
   |> Util.list_some
 
 
+let fci ?(frozen_core=true) mo_basis =
+  let mo_num = MOBasis.size mo_basis in
+  let ncore  = (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2 in
+  of_list (
+    if frozen_core then
+      List.concat [
+        Util.list_range 1 ncore
+        |> List.map (fun i -> Core i) ;
+        Util.list_range (ncore+1) mo_num
+        |> List.map (fun i -> Active i) 
+      ]
+    else
+      Util.list_range 1 mo_num
+      |> List.map (fun i -> Active i)
+  )
+

MOBasis/MOClass.mli

@@ -12,6 +12,11 @@ type t
 (** Creation *)
 val of_list : mo_class list -> t
 
+val fci : ?frozen_core:bool -> MOBasis.t -> t
+(** Creates the MO classes for FCI calculations : all [Active]. The 
+    [n] lowest MOs are [Core] if [frozen_core = true].
+*)
+
 
 val core_mos : t -> int list
 (** Returns a list containing the indices of the core MOs. *)

SCF/HartreeFock.ml

@@ -1,7 +1,7 @@
 open Util
 
 let make ?guess simulation =
-  if (Simulation.electrons simulation).Electrons.multiplicity = 1 then
+  if Electrons.multiplicity @@ Simulation.electrons simulation  = 1 then
     RHF.make ?guess simulation
   else
     invalid_arg "UHF or ROHF not implemented"

SCF/RHF.ml

@@ -12,7 +12,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
 
   (* Number of occupied MOs *)
   let nocc =
-    (Si.electrons simulation).El.n_alpha
+    El.n_alfa @@ Si.electrons simulation
   in
 
   let nuclear_repulsion = 

Utils/Electrons.ml

@@ -1,8 +1,8 @@
-(** Number of alpha and beta electrons *)
+(** Number of {% $\alpha$ %} and {% $\beta$ %} electrons *)
 
 type t = {
-  n_alpha : int ;
-  n_beta  : int ;
+  n_alfa : int ;
+  n_beta : int ;
   multiplicity : int;
 }
 
@@ -15,13 +15,16 @@ let make ?multiplicity:(multiplicity=1) ?charge:(charge=0) nuclei =
   let negative_charges = charge - positive_charges in
   let n_elec = - negative_charges in
   let n_beta = ((n_elec - multiplicity)+1)/2 in
-  let n_alpha = n_elec - n_beta in
-  let result = { n_alpha ; n_beta ; multiplicity } in
-  if multiplicity <> (n_alpha - n_beta)+1 then
+  let n_alfa = n_elec - n_beta in
+  let result = { n_alfa ; n_beta ; multiplicity } in
+  if multiplicity <> (n_alfa - n_beta)+1 then
     invalid_arg (__FILE__^": make");
   result
 
 let charge e =
-  - (e.n_alpha + e.n_beta)
+  - (e.n_alfa + e.n_beta)
   |> Charge.of_int
 
+let n_alfa t = t.n_alfa
+let n_beta t = t.n_beta
+let multiplicity t = t.multiplicity

Utils/Electrons.mli

@@ -1,11 +1,6 @@
 (** Information related to electrons. *)
 
-type t = {
-  n_alpha      : int ;  (** Number of alpha electrons *)
-  n_beta       : int ;  (** Number of beta  electrons *)
-  multiplicity : int ;  (** Spin multiplicity: {% $2S+1$ %} *)
-}
-
+type t 
 
 val make : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t
 (** Creates the data relative to electrons for a molecular system
@@ -19,3 +14,12 @@ val make : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t
 val charge : t -> Charge.t
 (** Sum of the charges of the electrons. *)
 
+val n_alfa : t -> int
+(** Number of alpha electrons *)
+
+val n_beta : t -> int
+(** Number of beta  electrons *)
+
+val multiplicity : t -> int 
+(** Spin multiplicity: {% $2S+1$ %} *)
+

Utils/Util.ml

@@ -298,8 +298,8 @@ let pp_matrix ppf m =
   aux 1 cols
 
 
-let pp_bitstring ppf bs =
-  String.init (Z.numbits bs) (fun i -> if (Z.testbit bs i) then '+' else '-')
+let pp_bitstring n ppf bs =
+  String.init n (fun i -> if (Z.testbit bs i) then '+' else '-')
     |> Format.fprintf ppf "@[<h>%s@]"
 
 

Utils/Util.mli

@@ -156,8 +156,8 @@ val pp_float_2darray : Format.formatter -> float array array -> unit
 ]}
 *)
 
-val pp_bitstring : Format.formatter -> Z.t -> unit
-(** Example: [ +++++------+-- ] *)
+val pp_bitstring : int -> Format.formatter -> Z.t -> unit
+(** Example: [ pp_bitstring 14 ppf z -> +++++------+-- ] *)
 
 val pp_matrix : Format.formatter -> Lacaml.D.mat -> unit