Improved arbitrary space

CI/CI.ml

@@ -1,6 +1,6 @@
 open Lacaml.D
 
-module Ds = Determinant_space
+module Ds = DeterminantSpace
 
 type t =
   {
@@ -53,24 +53,165 @@ let h_ij mo_basis ki kj =
 
 let create_matrix_arbitrary f det_space =
   lazy (
-    let det = 
+    let ndet = Ds.size det_space in
+    let data = 
       match Ds.determinants det_space with
-      | Ds.Arbitrary _ -> Ds.determinants_array det_space
+      | Ds.Arbitrary a -> a
       | _ -> assert false
     in
+    let det_alfa    = data.Ds.det_alfa
+    and det_beta    = data.Ds.det_beta
+    and det         = data.Ds.det
+    and index_start = data.Ds.index_start
+    in
 
-    let ndet = Ds.size det_space in
 
-    let v = Vec.make0 ndet in
+    (** Array of (list of singles, list of doubles) in the beta spin *)
+    let degree_bb = 
+      Array.map (fun det_i ->
+          let deg = Spindeterminant.degree det_i in
+          let doubles =
+            Array.mapi (fun i det_j -> 
+                let d = deg det_j in
+                if d < 3 then
+                  Some (i,d,det_j)
+                else
+                  None
+              ) det_beta
+            |> Array.to_list
+            |> Util.list_some
+          in
+          let singles = 
+            List.filter (fun (i,d,det_j) -> d < 2) doubles
+            |> List.map (fun (i,_,det_j) -> (i,det_j))
+          in
+          let doubles = 
+            List.map (fun (i,_,det_j) -> (i,det_j)) doubles
+          in
+          (singles, doubles)
+        ) det_beta
+    in
 
-    Array.init ndet
-      (fun i -> let ki = det.(i) in
-        Printf.eprintf "%8d / %8d\r%!" i ndet;
-        let j = ref 1 in
-        Ds.determinant_stream det_space
-        |> Stream.iter (fun kj -> v.{!j} <- f ki kj ; incr j);
-        Vector.sparse_of_vec v)
-    |> Matrix.sparse_of_vector_array
+    let task (i,i_dets) = 
+        let shift = index_start.(i) in
+
+        let result =
+          Array.init (index_start.(i+1) - shift)
+            (fun _ -> [])
+        in
+
+        (** Update function when ki and kj are connected *)
+        let update i j ki kj =
+          let x = f ki kj in
+          if abs_float x > Constants.epsilon then
+            result.(i - shift) <- (j, x) :: result.(i - shift) ;
+        in
+
+        let i_alfa = det_alfa.(i) in
+        let deg_a = Spindeterminant.degree i_alfa in
+
+        Array.iteri (fun j j_dets ->
+            let j_alfa = det_alfa.(j) in
+            let degree_a = deg_a j_alfa in
+
+            begin
+              match degree_a with
+              | 2 -> 
+              Array.iteri (fun i' i_b ->
+                  try
+                    Array.iteri (fun j' j_b ->
+                      if j_b >= i_b then 
+                      ( if j_b = i_b then
+                        ( let i_beta = det_beta.(i_b) in
+                          let ki = Determinant.of_spindeterminants i_alfa i_beta in
+                          let kj = Determinant.of_spindeterminants j_alfa i_beta in
+                          update (index_start.(i) + i')
+                                (index_start.(j) + j' + 1) ki kj);
+                        raise Exit)
+                      ) j_dets
+                    with Exit -> ()
+                ) i_dets 
+              | 1 -> 
+              Array.iteri (fun i' i_b ->
+                  let i_beta = det_beta.(i_b) in
+                  let ki = Determinant.of_spindeterminants i_alfa i_beta in
+                  let singles, _ = degree_bb.(i_b) in 
+                  let rec aux singles j' =
+                    match singles with
+                    | [] -> ()
+                    | (js, j_beta)::r_singles ->
+                      begin
+                        match compare js j_dets.(j') with
+                        | -1 -> aux r_singles j'
+                        |  0 ->
+                            let kj =
+                              Determinant.of_spindeterminants j_alfa j_beta
+                            in (update
+                                (index_start.(i) + i') (index_start.(j) + j' + 1)
+                                ki kj;
+                                aux r_singles (j'+1);)
+                        |  1 -> if (j' < Array.length j_dets) then aux singles (j'+1)
+                        |  _ -> assert false
+                      end
+                  in aux singles 0 
+                ) i_dets 
+              | 0 -> 
+              Array.iteri (fun i' i_b ->
+                  let i_beta = det_beta.(i_b) in
+                  let ki = Determinant.of_spindeterminants i_alfa i_beta in
+                  let _, doubles = degree_bb.(i_b) in 
+                  let rec aux doubles j' =
+                    match doubles with
+                    | [] -> ()
+                    | (js, j_beta)::r_doubles ->
+                      begin
+                        match compare js j_dets.(j') with
+                        | -1 -> aux r_doubles j'
+                        |  0 ->
+                            let kj =
+                              Determinant.of_spindeterminants j_alfa j_beta
+                            in (update
+                                (index_start.(i) + i') (index_start.(j) + j' + 1)
+                                ki kj;
+                                aux r_doubles (j'+1);)
+                        |  1 -> if (j' < Array.length j_dets) then aux doubles (j'+1)
+                        |  _ -> assert false
+                      end
+                  in aux doubles 0 
+                ) i_dets 
+              | _ -> ();
+            end
+          ) det;
+          let r =
+            Array.map (fun l -> 
+              List.rev l
+              |> Vector.sparse_of_assoc_list ndet
+            ) result
+          in (i,r)
+      in
+
+
+    let result =
+      if Parallel.master then
+        Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet [])
+      else
+        Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet [])
+    in
+
+    let n_det_alfa = 
+      Array.length det_alfa
+    in
+    Array.mapi (fun i i_dets -> i, i_dets) det
+    |> Array.to_list
+    |> Stream.of_list
+    |> Farm.run ~ordered:false ~f:task
+    |> Stream.iter (fun (k, r) ->
+        let shift = index_start.(k) in
+        let det_k = det.(k) in
+        Array.iteri (fun j r_j -> result.(shift+det_k.(j)) <- r_j) r;
+        Printf.eprintf "%8d / %8d\r%!" (k+1) n_det_alfa;
+      ) ;
+    Matrix.sparse_of_vector_array result
   )
 
 
@@ -117,13 +258,17 @@ let create_matrix_spin f det_space =
     in
 
     let task (i,i_alfa) = 
-        let result = Array.init n_beta (fun _ -> []) in
+        let result =
+          Array.init n_beta (fun _ -> [])
+        in
+
         (** Update function when ki and kj are connected *)
         let update i j ki kj =
           let x = f ki kj in
           if abs_float x > Constants.epsilon then
             result.(i) <- (j, x) :: result.(i) ;
         in
+
         let j = ref 1 in
         let deg_a = Spindeterminant.degree i_alfa in
         List.iter (fun j_alfa ->
@@ -198,7 +343,7 @@ let make ?(n_states=1) det_space =
     let mo_basis = Ds.mo_basis det_space in
 
     (* While in a sequential region, initiate the parallel
-       4-idx transformation
+       4-idx transformation to avoid nested parallel jobs
        *)
     ignore @@ MOBasis.two_e_ints mo_basis;
 

CI/DeterminantSpace.ml

@@ -32,7 +32,7 @@ type t =
   }
 
 
-module Ss = Spindeterminant_space
+module Ss = SpindeterminantSpace
 
 let n_alfa   t = t.n_alfa
 let n_beta   t = t.n_beta
@@ -45,7 +45,7 @@ let size t =
   | Spin (a,b)  ->  (Array.length a) * (Array.length b)
   | Arbitrary a ->  
       let ndet_a = Array.length a.det_alfa in
-      a.index_start.(ndet_a - 1) + Array.length a.det.(ndet_a - 1)
+      a.index_start.(ndet_a)
  
 
 let determinant_stream t =
@@ -150,15 +150,27 @@ let fci_of_mo_basis ?(frozen_core=true) mo_basis =
     let det_alfa = Ss.spin_determinants det_a 
     and det_beta = Ss.spin_determinants det_b
     in
-    Spin (det_alfa, det_beta)
-    (*
     let n_det_beta = Array.length det_beta in
-    Arbitrary {
-      det_alfa ; det_beta ;
-      det = Array.make (Array.length det_alfa) (Array.init (Array.length det_beta) (fun i -> i));
-      index_start = Array.mapi (fun i _ -> i*n_det_beta) det_alfa;
-    }
+    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 }
 

Makefile.include

@@ -49,7 +49,7 @@ doc: QCaml.odocl
 	$(OCAMLBUILD) -ocamlc ocamlcp $*.byte -use-ocamlfind $(PKGS)
 
 clean: 
-	$(OCAMLBUILD) -clean # rm -rf _build $(ALL_EXE) *.native *.byte
+	rm QCaml.odocl && $(OCAMLBUILD) -clean
 
 debug: run_integrals.native
 	./debug.sh

run_fci.ml

@@ -56,13 +56,13 @@ let () =
   in
 
   let space =
-    Determinant_space.fci_of_mo_basis ~frozen_core:false mos
+    DeterminantSpace.fci_of_mo_basis ~frozen_core:false mos
   in
   let ci = CI.make space in
   Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
   (*
   let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in
-  Util.list_range 1 (Determinant_space.size space)
+  Util.list_range 1 (DeterminantSpace.size space)
   |> List.iter (fun i -> Format.printf "@[%f@]@;" s2.{i,i});
   *)