diff --git a/CI/CI.ml b/CI/CI.ml index 3a97d7c..cad76d5 100644 --- a/CI/CI.ml +++ b/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 diff --git a/CI/DeterminantSpace.ml b/CI/DeterminantSpace.ml index e9a195b..d9cb913 100644 --- a/CI/DeterminantSpace.ml +++ b/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 } diff --git a/MOBasis/MOClass.ml b/MOBasis/MOClass.ml index b22c17b..6959aab 100644 --- a/MOBasis/MOClass.ml +++ b/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 diff --git a/MOBasis/MOClass.mli b/MOBasis/MOClass.mli index 6abd6f2..adf5392 100644 --- a/MOBasis/MOClass.mli +++ b/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 diff --git a/Utils/Util.ml b/Utils/Util.ml index afe94fa..7672584 100644 --- a/Utils/Util.ml +++ b/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} *) diff --git a/Utils/Util.mli b/Utils/Util.mli index 6f275e2..fba2036 100644 --- a/Utils/Util.mli +++ b/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