Accelerated Haa

CI/CI.ml

@@ -1,4 +1,5 @@
 open Lacaml.D
+open Util
 
 module Ds = DeterminantSpace
 module Sd = Spindeterminant
@@ -418,6 +419,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
         | MOClass.Core     _ -> false
         | _ -> true
       )
+    |> Array.of_list
   in
 
   let psi0 =
@@ -429,7 +431,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     Array.init (Ds.size det_space) (fun i ->
         Stream.next stream, psi0.{i+1,1})
   in
-    (*
+
   let is_internal =
     let m l =
       List.fold_left (fun accu i ->
@@ -455,20 +457,10 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
         in
         Z.logand neg_active_mask beta = occ_mask
   in
-        *)
 
 
   let symmetric = i_o1_alfa == alfa_o2_i 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 det_contribution i =
 
@@ -487,7 +479,9 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     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
+        | j ->
+          if Determinant.degree (fst psi0.(i)) (fst psi0.(j)) < 4
+          then
             aux (j::accu) (j-1)
           else
             aux accu (j-1)
@@ -522,13 +516,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     let det_i = fst psi0.(i) in
     let w_alfa = w_alfa det_i in
 
-
     let same_spin = 
       List.fold_left (fun accu spin ->
           accu +.
-          List.fold_left (fun accu particle ->
+          Array.fold_left (fun accu particle ->
               accu +.
-              List.fold_left (fun accu hole ->
+              Array.fold_left (fun accu hole ->
                   if hole = particle then accu else
                     let alfa =
                       Determinant.single_excitation spin hole particle det_i
@@ -541,12 +534,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                       in
 
                       let double =
-                        List.fold_left (fun accu particle' ->
+                        Array.fold_left (fun accu particle' ->
                             if particle' > particle || particle' = hole then
                               accu
                             else
                               accu +.
-                              List.fold_left (fun accu hole' ->
+                              Array.fold_left (fun accu hole' ->
                                   if hole' = particle' || hole' = particle || hole' < hole then
                                     accu
                                   else
@@ -570,9 +563,9 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     in
 
     let opposite_spin =
-      List.fold_left (fun accu particle ->
+      Array.fold_left (fun accu particle ->
           accu +.
-          List.fold_left (fun accu hole ->
+          Array.fold_left (fun accu hole ->
               if hole = particle then accu else
                 let alfa =
                   Determinant.single_excitation Spin.Alfa hole particle det_i
@@ -580,9 +573,9 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                 if Determinant.is_none alfa then accu else
 
                   let double =
-                    List.fold_left (fun accu particle' ->
+                    Array.fold_left (fun accu particle' ->
                         accu +.
-                        List.fold_left (fun accu hole' ->
+                        Array.fold_left (fun accu hole' ->
                             if hole' = particle' then accu else
                               let alfa = 
                                 Determinant.double_excitation
@@ -616,9 +609,48 @@ let pt2_en ci =
 
   let i_o1_alfa = h_ij mo_basis in
 
-  let w_alfa _ =
+  let w_alfa det_i =
+    let one_e, two_e = h_integrals mo_basis in
+    let fock_diag_alfa, fock_diag_beta = 
+      Ds.fock_diag ci.det_space det_i
+    in
+    let h_aa alfa =
+      match Excitation.of_det det_i alfa with
+      | Excitation.Double (_,
+            {hole = h ; particle = p ; spin = s }, 
+            {hole = h'; particle = p'; spin = s'}) -> 
+              let fock_diag1 =
+                if s = Spin.Alfa then
+                  fock_diag_alfa
+                else
+                  fock_diag_beta
+              in
+              let fock_diag2 =
+                if s' = Spin.Alfa then
+                  fock_diag_alfa
+                else
+                  fock_diag_beta
+              in
+              fock_diag1.(0) -. fock_diag1.(h)
+              +. (fock_diag1.(p ) -. two_e h p h p s s)
+              -. (fock_diag2.(h') -. two_e h h' h h' s s' +. two_e p h' p h' s s')
+              +. (fock_diag2.(p') -. two_e h p' h p' s s'
+                  +. two_e p p' p p' s s' -. two_e h' p' h' p' s' s')
+      | Excitation.Single (_,
+            {hole = h ; particle = p ; spin = s }) -> 
+              let fock_diag =
+                if s = Spin.Alfa then
+                  fock_diag_alfa
+                else
+                  fock_diag_beta
+              in
+              fock_diag.(0) -. fock_diag.(h)
+              +. (fock_diag.(p) -. two_e h p h p s s)
+              |> ignore;
+              h_ij mo_basis alfa alfa
+      | _ -> assert false
+    in
     let e0 = e0.{1} in 
-    let h_aa alfa = h_ij mo_basis alfa alfa in
     fun alfa ->
       1. /. (e0 -. h_aa alfa)
   in

CI/Determinant.mli

@@ -58,7 +58,6 @@ val degree : t -> t -> int
 (** Returns degree of excitation between two determinants. *)
 
 
-
 (** {1 Creators} *)
 
 val of_spindeterminants : Spindeterminant.t -> Spindeterminant.t -> t

CI/DeterminantSpace.ml

@@ -154,6 +154,78 @@ let determinant t i =
     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 = MOBasis.size mo_basis in
+  let one_e_ints  = MOBasis.one_e_ints mo_basis
+  and two_e_ints  = MOBasis.two_e_ints mo_basis
+  in
+  let two_e i j k l = ERI.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 cas_of_mo_basis mo_basis n m =
   let s = MOBasis.simulation mo_basis in
   let e = Simulation.electrons s in

CI/DeterminantSpace.mli

@@ -44,6 +44,11 @@ val determinant_stream : t -> Determinant.t Stream.t
 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 : ?frozen_core:bool -> MOBasis.t -> t
 (** Creates a space of all possible ways to put [n_alfa] electrons in the {% $\alpha$ %}

CI/Spindeterminant.ml

@@ -113,7 +113,6 @@ let of_list l =
   |> List.fold_left (fun accu p -> creation p accu) vac
 
 
-
 let rec to_list = function
   | None -> []
   | Some spindet ->
@@ -125,6 +124,11 @@ let rec to_list = function
     in aux [] spindet.bitstring
 
 
+let n_electrons = function
+  | Some t -> Z.popcount t.bitstring
+  | None -> 0
+
+
 let pp_spindet n ppf = function
   | None -> Format.fprintf ppf "@[<h>None@]"
   | Some s -> 

CI/Spindeterminant.mli

@@ -57,6 +57,10 @@ 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 -> Z.t -> t