diff --git a/CI/CI.ml b/CI/CI.ml index e222865..91c4e45 100644 --- a/CI/CI.ml +++ b/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 = lazy ( - let ntasks = int_of_float @@ sqrt @@ float_of_int ndet in - List.init ntasks (fun i -> - let m = - if i = (ntasks-1) then - (ndet - i*ntasks) - else - ntasks - in - List.init m (fun j -> i*ntasks + j) - ) - |> 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 + let m_H = + let m_H_arbitrary det = lazy ( + let v = Vec.make0 ndet 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} <- h_ij mo_basis ki kj ; incr j); + Vector.sparse_of_vec v) + |> Matrix.sparse_of_vector_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 ( - let v = Vec.make0 ndet 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} <- h_ij mo_basis ki kj ; incr j); - Vector.sparse_of_vec v) - |> Matrix.sparse_of_vector_array - ) + match Ds.determinants det_space with + | Ds.Arbitrary a -> m_H_arbitrary a + | Ds.Spin (a,b) -> m_H_spin a b in + let m_S2 = lazy ( - let v = Vec.make0 ndet in - Array.init ndet (fun i -> let ki = det.(i) in - Array.iteri (fun j kj -> - v.{j+1} <- CIMatrixElement.make_s2 ki kj) det; - Vector.sparse_of_vec v) - |> Matrix.sparse_of_vector_array + match Ds.determinants det_space with + | Ds.Spin (a,b) -> failwith "Not implemented" + | Ds.Arbitrary det -> + let v = Vec.make0 ndet in + Array.init ndet (fun i -> let ki = det.(i) in + 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 - (* - Parallel.broadcast @@ - lazy (Util.diagonalize_symm m_H) - *) + let m_H = + Lazy.force m_H + in 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 } diff --git a/CI/Determinant_space.ml b/CI/Determinant_space.ml index 3eefe87..6fcee7c 100644 --- a/CI/Determinant_space.ml +++ b/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 -> - if i < imax then Some a.(i) else None) + Stream.from (fun i -> + 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) -> diff --git a/CI/Determinant_space.mli b/CI/Determinant_space.mli index c12b456..bd7b42a 100644 --- a/CI/Determinant_space.mli +++ 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