Parallelized PT2

CI/CI.ml

@@ -24,6 +24,8 @@ let m_S2         t = Lazy.force t.m_S2
 
 let eigensystem  t = Lazy.force t.eigensystem 
 
+let mo_class     t = DeterminantSpace.mo_class t.det_space
+
 let eigenvectors t = 
   let (x,_) = eigensystem t in x
 
@@ -394,7 +396,9 @@ let make ?(n_states=1) det_space =
       Matrix.mm ~transa:`T m_H psi
     in
     let eigenvectors, eigenvalues = 
+      Parallel.broadcast (lazy (
       Davidson.make ~threshold:1.e-6  ~n_states diagonal matrix_prod
+      ))
     in
     let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in
     eigenvectors, eigenvalues
@@ -406,27 +410,14 @@ let make ?(n_states=1) det_space =
 
 
 let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
-    i_o1_alfa  alfa_o2_i  w_alfa =
+    list_holes  list_particles i_o1_alfa  alfa_o2_i  w_alfa =
 
-  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 mo_class in
-    Util.list_range 1 (MOBasis.size mo_basis)
-    |> List.filter (fun i -> match cls.(i) with
-        | MOClass.Inactive _
-        | MOClass.Active _
-        | MOClass.Virtual _ -> true
-        | MOClass.Auxiliary _
-        | MOClass.Deleted _
-        | MOClass.Core     _ -> false
-      )
-    |> Array.of_list
+  let list_holes     = Array.of_list list_holes 
+  and list_particles = Array.of_list list_particles
   in
 
   let psi0 =
-    let psi0, _ = Lazy.force eigensystem in
+    let psi0, _ = Parallel.broadcast eigensystem in
 
     let stream = 
       Ds.determinant_stream det_space
@@ -441,6 +432,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
           let j = i-1 in Z.(logor accu (shift_left one j))
         ) Z.zero l
     in
+    let mo_class      = DeterminantSpace.mo_class det_space 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
@@ -556,12 +548,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                                         accu
                                       else
                                         accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa 
-                                ) 0. mo_indices
-                          ) 0. mo_indices
+                                ) 0. list_holes
+                          ) 0. list_particles
                       in
                       accu +. single +. double
-                ) 0. mo_indices
-            ) 0. mo_indices
+                ) 0. list_holes
+            ) 0. list_particles
         ) 0. [ Spin.Alfa ; Spin.Beta ]
     in
 
@@ -590,25 +582,28 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                                 accu
                               else
                                 accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa 
-                          ) 0. mo_indices
-                      ) 0. mo_indices
+                          ) 0. list_holes
+                      ) 0. list_particles
                   in
                   accu +. double
-            ) 0. mo_indices
-        ) 0. mo_indices
+            ) 0. list_holes
+        ) 0. list_particles
     in
     same_spin +. opposite_spin
   in
 
-  Array.mapi (fun i (_,c_i) -> det_contribution i) psi0
-  |> Array.fold_left (+.) 0.
-
+  let result = 
+    Util.stream_range 0 (Array.length psi0 - 1)
+    |> Farm.run ~ordered:false ~f:det_contribution
+    |> Util.stream_fold (+.) 0.
+  in
+  Parallel.broadcast (lazy result)
 
 
 let pt2_en ci =
 
   let mo_basis  = Ds.mo_basis ci.det_space in
-  let _psi0, e0 = Lazy.force ci.eigensystem in
+  let _psi0, e0 = Parallel.broadcast ci.eigensystem in
 
   let i_o1_alfa = h_ij mo_basis in
 
@@ -658,7 +653,16 @@ let pt2_en ci =
       1. /. (e0 -. h_aa alfa)
   in
 
-  second_order_sum  ci  i_o1_alfa  i_o1_alfa  w_alfa
+  let mo_class = mo_class ci in
+  let list_holes = List.concat
+    [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
+  and list_particles = List.concat 
+    [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
+  in
+
+  second_order_sum  ci  list_holes  list_particles
+    i_o1_alfa  i_o1_alfa  w_alfa
+
 
 
 
@@ -679,7 +683,15 @@ let pt2_mp ci =
     | _ -> assert false
   in
 
-  second_order_sum  ci  i_o1_alfa  i_o1_alfa  w_alfa
+  let mo_class = mo_class ci in
+  let list_holes = List.concat
+    [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
+  and list_particles = List.concat 
+    [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
+  in
+
+  second_order_sum  ci  list_holes  list_particles
+    i_o1_alfa  i_o1_alfa  w_alfa
 
 
 let variance ci =
@@ -690,6 +702,14 @@ let variance ci =
 
   let w_alfa _ _ = 1. in
 
-  second_order_sum  ci  i_o1_alfa  i_o1_alfa  w_alfa
+  let mo_class = mo_class ci in
+  let list_holes = List.concat
+    [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
+  and list_particles = List.concat 
+    [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
+  in
+
+  second_order_sum  ci  list_holes  list_particles
+    i_o1_alfa  i_o1_alfa  w_alfa
 
 

CI/DeterminantSpace.ml

@@ -245,8 +245,9 @@ let spin_of_mo_basis mo_basis f =
     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;
+    if Parallel.master then
+      Format.printf "Number of determinants : %d %d %d\n%!" 
+        n_det_alfa n_det_beta ndet;
     Spin (det_alfa, det_beta)
   in
   { n_alfa ; n_beta ; mo_class ; mo_basis ; determinants }
@@ -270,8 +271,9 @@ let arbitrary_of_mo_basis mo_basis f =
     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;
+    if Parallel.master then
+      Format.printf "Number of determinants : %d %d %d\n%!" 
+        n_det_alfa n_det_beta ndet;
 
     let det = Array.make n_det_alfa
         (Array.init n_det_beta (fun i -> i))
@@ -279,8 +281,9 @@ let arbitrary_of_mo_basis mo_basis f =
     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;
+    if Parallel.master then
+      Format.printf "Number of determinants : %d %d %d\n%!" 
+        n_det_alfa n_det_beta ndet;
     Arbitrary {
       det_alfa ; det_beta ; det ; index_start
     }

MOBasis/MOClass.ml

@@ -64,7 +64,7 @@ let deleted_mos t =
   |> Util.list_some
 
 
-let virtual_mos t =
+let auxiliary_mos t =
   List.map (fun x ->
       match x with 
       | Auxiliary i -> Some i

MOBasis/MOClass.mli

@@ -42,6 +42,9 @@ val inactive_mos : t -> int list
 val deleted_mos : t -> int list
 (** Returns a list containing the indices of the deleted MOs. *)
 
+val auxiliary_mos : t -> int list
+(** Returns a list containing the indices of the auxiliary MOs. *)
+
 val mo_class_array : t -> mo_class array
 (** Returns an array [a] such that [a.(i)] returns the class of MO [i].
     As the MO indices start from [1], the array has an extra zero entry

Utils/Util.ml

@@ -239,6 +239,16 @@ let stream_to_list stream =
   in aux []
 
 
+let stream_fold f init stream =
+  let rec aux accu =
+    try
+      let element = Stream.next stream in
+      let new_accu = f accu element in
+      aux new_accu
+    with Stream.Failure -> accu
+  in
+  aux init
+
 
 (** {2 Linear algebra} *)
 

Utils/Util.mli

@@ -104,6 +104,10 @@ val stream_range : int -> int -> int Stream.t
 val stream_to_list : 'a Stream.t -> 'a list
 (** Read a stream and put items in a list. *)
 
+val stream_fold : ('a -> 'b -> 'a) -> 'a -> 'b Stream.t -> 'a
+(** Apply a fold to the elements of the stream. *)
+
+
 (** {2 Linear algebra } *)
 
 val normalize : Vec.t -> Vec.t