PT2 is not working yet

CI/CI.ml

@@ -385,12 +385,14 @@ let make ?(n_states=1) det_space =
 
 
 
-(*
 let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } =
 
+  let mo_basis = Ds.mo_basis det_space in
+    
+  let mo_class = DeterminantSpace.mo_class det_space in
   let mo_indices =
-    let cls = MOClass.mo_class_array (DeterminantSpace.mo_class det_space) in
+    let cls = MOClass.mo_class_array mo_class in
-    Util.list_range 1 (Ds.mo_basis det_space |> MOBasis.size)
+    Util.list_range 1 (MOBasis.size mo_basis)
     |> List.filter (fun i -> match cls.(i) with
       | MOClass.Deleted  _ 
       | MOClass.Core     _ -> false
@@ -398,6 +400,7 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } =
       )
   in
 
+(* Only the gournd state is computed here *)
   let psi0, e0 = Lazy.force eigensystem in
 
   let psi0 =
@@ -405,44 +408,128 @@ let pt2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states } =
       Ds.determinant_stream det_space
     in
     Array.init (Ds.size det_space) (fun i ->
-      (Stream.next stream, psi0.{i,1})
+      Stream.next stream, psi0.{i+1,1})
   in
   let e0 = e0.{1} in 
 
-  let det_contribution i =
+    (*
-    let psi_filtered =
+  let is_internal =
-      List.filter (fun (det, _) ->
+    let m l =
-          Determinant.degree det i < 4) psi0
+      List.fold_left (fun accu i ->
+          let j = i-1 in Z.(logor accu (shift_left one j))
+        ) Z.zero l
     in
-    List.fold_left (fun accu spin1 ->
+    let active_mask   = m (MOClass.active_mos   mo_class) in
-        accu +.
+    let occ_mask      = m (MOClass.core_mos     mo_class) in
-        List.fold_left (fun accu particle ->
+    let inactive_mask = m (MOClass.inactive_mos mo_class) in
-            accu +.
+    let occ_mask = Z.logor occ_mask inactive_mask in
-            List.fold_left (fun accu hole ->
+    let neg_active_mask = Z.lognot active_mask in
-                let alfa =
+    fun a ->
-                  Determinant.single_excitation spin1 hole particle i
+      let alfa =
-                in
+        Determinant.alfa a
-                if Determinant.is_none alfa then
+        |> Spindeterminant.bitstring
-                  accu
+      in
-                else
+      if Z.logand neg_active_mask alfa <> occ_mask then
-                  let psi_h_alfa =
+        false
-                    Array.fold_left (fun (det, coef) ->
+      else
-                        coef *. (h_ij det alfa)
+        let beta =
-                      ) 0. psi_filtered
+          Determinant.beta a
-                  in
+          |> Spindeterminant.bitstring
-                  let h_aa = h_ij alfa alfa in
+        in
-                  accu +. psi_h_alfa *. psi_h_alfa / (e0 - h_aa)
+        Z.logand neg_active_mask beta = occ_mask
-              ) 0. mo_indices
-          ) 0. mo_indices
-      ) 0. [ Spin.Alfa ; Spin.Beta ]
-
   in
+        *)
+
+    
+  let det_contribution i =
+
+    let psi_filtered_idx =
+      let rec aux accu =  function
+      | j when j <= i -> List.rev accu
+      | j -> if Determinant.degree (fst psi0.(i)) (fst psi0.(j)) < 4 then
+                aux (j::accu) (j-1)
+             else
+                aux accu (j-1)
+      in aux [] (Array.length psi0 - 1)
+    in
+
+    let psi_filtered = 
+      List.map (fun i -> psi0.(i)) psi_filtered_idx
+    in
+
+    let psi_h alfa =
+      let hij = h_ij mo_basis alfa in
+      List.fold_left (fun accu (det, coef) ->
+          accu +. coef *. (hij det)) 0. psi_filtered
+    in
+
+    let is_internal alfa =
+        let rec aux = function
+        | -1 -> false
+        | j  -> Determinant.degree (fst psi0.(j)) alfa = 0 || aux (j-1)
+        in
+        aux (Array.length psi0 - 1)
+    in
+
+    let already_generated alfa =
+      if is_internal alfa then
+        true
+      else
+        let rec aux = function
+        | -1 -> false
+        | j  -> Determinant.degree (fst psi0.(j)) alfa <= 2 || aux (j-1)
+        in
+        aux (i-1)
+    in
+
+    let det_i = fst psi0.(i) in
+
+    List.fold_left (fun accu particle ->
+        accu +.
+        List.fold_left (fun accu hole ->
+            if hole = particle then accu else
+              accu +.
+              List.fold_left (fun accu spin ->
+                  let alfa =
+                    Determinant.single_excitation spin hole particle det_i
+                  in
+                  if Determinant.is_none alfa then accu else
+
+                    let single = 
+                      if already_generated alfa then 0. else
+                        let h_aa = h_ij mo_basis alfa alfa in
+                        let psi_h_alfa = psi_h alfa in
+                        psi_h_alfa *. psi_h_alfa /. (e0 -. h_aa)
+                    in
+
+                    let double =
+                      List.fold_left (fun accu particle' ->
+                          accu +.
+                          List.fold_left (fun accu hole' ->
+                              accu +.
+                              List.fold_left (fun accu spin' ->
+                                  let alfa = 
+                                      Determinant.double_excitation
+                                        spin  hole  particle
+                                        spin' hole' particle' det_i
+                                  in
+                                  if Determinant.is_none alfa ||
+                                    already_generated alfa then
+                                    accu
+                                  else
+                                      let h_aa = h_ij mo_basis alfa alfa in
+                                      let psi_h_alfa = psi_h alfa in
+                                      accu +.  psi_h_alfa *. psi_h_alfa /. (e0 -. h_aa)
+                                  ) 0. [ Spin.Alfa ; Spin.Beta ]
+                              ) 0. mo_indices
+                          ) 0. mo_indices
+                    in
+                    accu +. single +. double
+                ) 0. [ Spin.Alfa ; Spin.Beta ]
+          ) 0. mo_indices
+      ) 0. mo_indices
+    in
+
+  Array.mapi (fun i (_,c_i) -> c_i *. det_contribution i) psi0
+  |> Array.fold_left (+.) 0.
 
-  match det_space with
-  | Ds.Arbitrary -> assert false
-  | Ds.Spin alfa_dets beta_dets ->
-    Array.map ( fun alfa ->
-        Array.map ( fun beta ->
-          ) beta_dets
-      ) alfa_dets
- *) 

CI/DeterminantSpace.ml

@@ -154,6 +154,46 @@ let determinant t i =
     Determinant.of_spindeterminants alfa beta
 
 
+let cas_of_mo_basis mo_basis n m =
+  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.cas_of_mo_basis mo_basis n_alfa n m
+  and det_b =
+    Ss.cas_of_mo_basis mo_basis n_beta n m
+  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)
+
+    (*
+    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 fci_of_mo_basis ?(frozen_core=true) mo_basis =
   let s = MOBasis.simulation mo_basis in
   let e = Simulation.electrons s in

CI/DeterminantSpace.mli

@@ -51,6 +51,10 @@ val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> t
     All other MOs are untouched.
 *)
 
+val cas_of_mo_basis : MOBasis.t -> int -> int -> t
+(** Creates a CAS(n,m) space of determinants.
+*)
+
 (** {2 Printing} *)
 
 val pp_det_space : Format.formatter -> t -> unit

CI/SpindeterminantSpace.ml

@@ -24,6 +24,8 @@ let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num =
   and active_mask = m (MOClass.active_mos  mo_class)
   in
   let neg_active_mask = Z.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 = 
     Util.bit_permtutations elec_num mo_num 
     |> List.filter (fun b -> Z.logand neg_active_mask b = occ_mask)
@@ -33,6 +35,30 @@ let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num =
   { elec_num ; mo_basis ; mo_class ; spin_determinants }
 
 
+let cas_of_mo_basis mo_basis elec_num n m =
+  let mo_num = MOBasis.size mo_basis in
+  let mo_class = MOClass.cas_sd mo_basis n m 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
+  in
+  let active_mask   = m (MOClass.active_mos   mo_class) in
+  let occ_mask      = m (MOClass.core_mos     mo_class) in
+  let inactive_mask = m (MOClass.inactive_mos mo_class) in
+  let occ_mask = Z.logor occ_mask inactive_mask in
+  let neg_active_mask = Z.lognot active_mask in
+  (* Here we generate the FCI space and filter out all the unwanted determinants.
+     This should be improved. *)
+  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 ; spin_determinants }
+
+
+
 let pp_spindet_space ppf t =
   Format.fprintf ppf "@[<v 2>[ ";
   Array.iteri (fun i d -> Format.fprintf ppf "@[<v>@[%8d@] @[%a@]@]@;" i

CI/SpindeterminantSpace.mli

@@ -25,8 +25,14 @@ val mo_basis : t -> MOBasis.t
 (** {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.
+(** Create a space of all possible ways to put [n_elec-ncore] electrons in the
-    All other MOs are untouched.
+    [Active] MOs.  All other MOs are untouched.
+*)
+
+val cas_of_mo_basis : MOBasis.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].
 *)