spin fci

CI/CI.ml

@@ -8,10 +8,13 @@ type t =
     m_H         : Matrix.t lazy_t ;
     m_S2        : Matrix.t lazy_t ;
     eigensystem : (Mat.t * Vec.t) lazy_t;
+    n_states    : int;
   }
 
 let det_space    t = t.det_space
 
+let n_states     t = t.n_states
+
 let m_H          t = Lazy.force t.m_H
 
 let m_S2         t = Lazy.force t.m_S2
@@ -49,108 +52,82 @@ let h_ij mo_basis ki kj =
 
 
 
-let make det_space = 
+let make ?(n_states=1) det_space = 
   let ndet = Ds.size det_space in
-  let det  = Ds.determinants det_space in
   let mo_basis = Ds.mo_basis det_space in
-  
 
-  (*
+  let m_H = 
-  let m_H = lazy (
+    let m_H_arbitrary det = lazy (
-    let ntasks = int_of_float @@ sqrt @@ float_of_int ndet in
+      let v = Vec.make0 ndet in
-    List.init ntasks (fun i -> 
+      Array.init ndet (fun i -> let ki = det.(i) in
-      let m = 
+        Printf.eprintf "%8d / %8d\r%!" i ndet;
-        if i = (ntasks-1) then
+        let j = ref 1 in
-          (ndet - i*ntasks) 
+        Ds.determinant_stream det_space
-        else
+        |> Stream.iter (fun kj ->
-          ntasks
+          v.{!j} <- h_ij mo_basis ki kj ; incr j);
-      in
+        Vector.sparse_of_vec v)
-      List.init m (fun j -> i*ntasks + j)
+      |> Matrix.sparse_of_vector_array)
-      )
-    |> Stream.of_list
-    |> Farm.run ~ordered:true
-      ~f:(fun l ->
-          Printf.eprintf "Start\n%!";
-      List.map (fun i ->
-          let ki = det.(i) in
-          Printf.eprintf "%d / %d\n%!" i ndet;
-          Array.init ndet (fun j -> let kj = det.(j) in
-                            h_ij mo_basis ki kj)
-          ) l)
-    |> Util.stream_to_list
-    |> List.concat
-    |> List.map Vec.of_array
-    |> Mat.of_col_vecs_list
-  ) in
-  *)
-
-  (* Parallel
-  let m_H = lazy (
-    let h =
-      if Parallel.master then
-        Array.make_matrix ndet ndet 0.
-        |> Mat.of_array
-      else
-        Array.make_matrix 1 1 0.
-        |> Mat.of_array
     in
-    List.init ndet (fun i -> i)
-    |> Stream.of_list
-    |> Farm.run ~ordered:false
-      ~f:(fun i ->
-          let ki = det.(i) in
-          Printf.eprintf "%8d / %8d\r%!" i ndet;
-          List.init ndet (fun j ->
-            let kj = det.(j) in
-            let x = h_ij mo_basis ki kj in
-            if x <> 0. then Some (i,j,x) else None
-          )
-          |> Util.list_some
-          )
-    |> Util.stream_to_list
-    |> List.iter (fun l -> if Parallel.master then
-      List.iter (fun (i,j,x) -> h.{i+1,j+1} <- x) l);
-    Parallel.broadcast (lazy h)
-  ) in
-  *)
 
+    let m_H_spin a b = lazy (
+      let v = Vec.make0 ndet in
+      let i = ref 0 in
+      let result = Array.make ndet (Vector.sparse_of_vec @@ Vec.make0 1) in
+      Array.iteri (fun ia i_alfa ->
+          Array.iteri (fun ib i_beta ->
+              Printf.eprintf "%8d / %8d\r%!" !i ndet;
+              let ki =
+                Determinant.of_spindeterminants i_alfa i_beta
+              in
+              let j = ref 0 in
+              Array.iteri (fun ja j_alfa ->
+                  Array.iteri (fun jb j_beta ->
+                      let kj =
+                        Determinant.of_spindeterminants j_alfa j_beta
+                      in
+                      incr j;
+                      v.{!j} <- h_ij mo_basis ki kj
+                    ) b;
+                ) a;
+              result.(!i) <- Vector.sparse_of_vec v;
+              incr i;
+            ) b
+        ) a;
+      Matrix.sparse_of_vector_array result
+    )
+    in
 
-  let m_H = lazy (
+    match Ds.determinants det_space with
-    let v = Vec.make0 ndet in
+    | Ds.Arbitrary a -> m_H_arbitrary a
-    Array.init ndet (fun i -> let ki = det.(i) in
+    | Ds.Spin  (a,b) -> m_H_spin a b 
-      Printf.eprintf "%8d / %8d\r%!" i ndet;
-      let j = ref 1 in
-      Ds.determinant_stream det_space
-      |> Stream.iter (fun kj ->
-        v.{!j} <- h_ij mo_basis ki kj ; incr j);
-      Vector.sparse_of_vec v)
-    |> Matrix.sparse_of_vector_array
-  )
   in
+
   let m_S2 = lazy (
-    let v = Vec.make0 ndet in
+    match Ds.determinants det_space with
-    Array.init ndet (fun i -> let ki = det.(i) in
+    | Ds.Spin  (a,b) -> failwith "Not implemented"
-      Array.iteri (fun j kj ->
+    | Ds.Arbitrary det -> 
-        v.{j+1} <- CIMatrixElement.make_s2 ki kj) det;
+      let v = Vec.make0 ndet in
-      Vector.sparse_of_vec v)
+      Array.init ndet (fun i -> let ki = det.(i) in
-    |> Matrix.sparse_of_vector_array
+        Array.iteri (fun j kj ->
+          v.{j+1} <- CIMatrixElement.make_s2 ki kj) det;
+        Vector.sparse_of_vec v)
+      |> Matrix.sparse_of_vector_array
   )
   in
+
   let eigensystem = lazy (
-    let m_H = Lazy.force m_H in
+    let m_H =
-    (*
+      Lazy.force m_H
-    Parallel.broadcast @@
+    in
-      lazy (Util.diagonalize_symm m_H)
-      *)
     let diagonal =
       Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i)
     in
     let matrix_prod psi = 
       Matrix.mm ~transa:`T m_H psi
     in
-    Davidson.make diagonal matrix_prod
+    Davidson.make ~n_states diagonal matrix_prod
   )
   in
-  { det_space ; m_H ; m_S2 ; eigensystem }
+  { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
 
 

CI/Determinant_space.ml

@@ -30,8 +30,8 @@ let determinant_stream t =
   let imax = size t in
   match t.determinants with
   | Arbitrary a -> 
-    Stream.from (fun i ->
+      Stream.from (fun i ->
-        if i < imax then Some a.(i) else None)
+          if i < imax then Some a.(i) else None)
   | Spin (a,b) -> 
     let na = Array.length a
     and nb = Array.length b in
@@ -49,7 +49,10 @@ let determinant_stream t =
           None)
 
 
-let determinants t =
+let determinants t = t.determinants
+
+
+let determinants_array t =
   match t.determinants with
   | Arbitrary a -> a
   | Spin (a,b)  -> 

CI/Determinant_space.mli

@@ -4,6 +4,10 @@ The determinant space in which we solve the Schrodinger equation.
 
 type t
 
+type determinant_storage =
+| Arbitrary of Determinant.t array
+| Spin      of (Spindeterminant.t array * Spindeterminant.t array)
+
 (** {1 Accessors} *)
 
 val n_alfa : t -> int
@@ -18,7 +22,10 @@ val mo_class : t -> MOClass.t
 val mo_basis : t -> MOBasis.t
 (** The MO basis on which the determinants are expanded. *)
 
-val determinants : t -> Determinant.t array
+val determinants : t -> determinant_storage
+(** All the determinants belonging to the space. *)
+
+val determinants_array : t -> Determinant.t array
 (** All the determinants belonging to the space. *)
 
 val determinant_stream : t -> Determinant.t Stream.t