From 8b49ac8f77d39c9efded49932676129b01af668d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 3 Oct 2019 16:58:15 +0200 Subject: [PATCH] Dressed integrals in F12 --- CI/F12CI.ml | 309 +++++---------------------------------- MOBasis/HF12.ml | 186 +++++++++++++++++++++++ Utils/FourIdxStorage.ml | 3 + Utils/FourIdxStorage.mli | 6 + run_fci_f12.ml | 1 + 5 files changed, 236 insertions(+), 269 deletions(-) create mode 100644 MOBasis/HF12.ml diff --git a/CI/F12CI.ml b/CI/F12CI.ml index ea87548..4363527 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -1,314 +1,83 @@ open Lacaml.D +module Ds = DeterminantSpace + type t = { mo_basis : MOBasis.t ; - aux_basis : MOBasis.t ; det_space : DeterminantSpace.t ; ci : CI.t ; + hf12_integrals : HF12.t ; eigensystem : (Mat.t * Vec.t) lazy_t; } let ci t = t.ci let mo_basis t = t.mo_basis let det_space t = t.det_space -let mo_class t = DeterminantSpace.mo_class @@ det_space t +let mo_class t = Ds.mo_class @@ det_space t let eigensystem t = Lazy.force t.eigensystem -let f12_integrals mo_basis = - let two_e_ints = MOBasis.f12_ints mo_basis in - ( (fun _ _ _ -> 0.), - (fun i j k l s s' -> - if (i=k && j<>l) || (j=l && i<>k) then - 0. - else - begin - let ijkl = F12.get_phys two_e_ints i j k l - in - (* - if s' = Spin.other s then - (* Minus sign because we swap spin variables - instead of orbital variables *) - 0.375 *. ijkl +. 0.125 *. ijlk - else - 0.25 *. (ijkl -. ijlk) - *) - if s' = Spin.other s then - ijkl - else - let ijlk = F12.get_phys two_e_ints i j l k - in - ijkl -. ijlk - end - ) ) - -let h_ij mo_basis ki kj = - let integrals = - List.map (fun f -> f mo_basis) - [ CI.h_integrals ] +let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = + let integrals = [ + let one_e _ _ _ = 0. in + let hf12_s, hf12_o = hf12_integrals in + let two_e i j k l s s' = + if s' = Spin.other s then + hf12_o.{i,j,k,l} + else + hf12_s.{i,j,k,l} + in + (one_e, two_e) + ] in - CIMatrixElement.make integrals ki kj - |> List.hd - - -let f_ij mo_basis ki kj = - let integrals = - List.map (fun f -> f mo_basis) - [ f12_integrals ] - in - CIMatrixElement.make integrals ki kj + CIMatrixElement.non_zero integrals deg_a deg_b ki kj |> List.hd -let hf_ij mo_basis ki kj = - let integrals = - List.map (fun f -> f mo_basis) - [ CI.h_integrals ; f12_integrals ] - in - CIMatrixElement.make integrals ki kj - - -let is_a_double det_space = - let mo_class = DeterminantSpace.mo_class det_space in - let mo_num = Array.length @@ MOClass.mo_class_array mo_class in - let m l = - List.fold_left (fun accu i -> - let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j) - ) (Bitstring.zero mo_num) l - in - let aux_mask = m (MOClass.auxiliary_mos mo_class) in - fun k -> - let alfa = - Determinant.alfa k - |> Spindeterminant.bitstring - in - let beta = - Determinant.beta k - |> Spindeterminant.bitstring - in - let a = Bitstring.logand aux_mask alfa - and b = Bitstring.logand aux_mask beta - in - match Bitstring.popcount a + Bitstring.popcount b with - | 2 -> true - | _ -> false - - -let p12 det_space = - let mo_class = DeterminantSpace.mo_class det_space in - let mo_num = Array.length @@ MOClass.mo_class_array mo_class in - let m l = - List.fold_left (fun accu i -> - let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j) - ) (Bitstring.zero mo_num) l - in - let aux_mask = m (MOClass.auxiliary_mos mo_class) in - let not_aux_mask = - Bitstring.(shift_left_one mo_num mo_num |> minus_one) - in - fun k -> - let alfa = - Determinant.alfa k - |> Spindeterminant.bitstring - in - let beta = - Determinant.beta k - |> Spindeterminant.bitstring - in - let a = Bitstring.logand aux_mask alfa - and b = Bitstring.logand aux_mask beta - in - match Bitstring.popcount a, Bitstring.popcount b with - | 2, 0 - | 0, 2 -> Some (Determinant.negate_phase k) - | 1, 1 -> Some (Determinant.of_spindeterminants - (Spindeterminant.of_bitstring @@ - Bitstring.(logor b (logand not_aux_mask alfa)) ) - (Spindeterminant.of_bitstring @@ - Bitstring.(logor a (logand not_aux_mask beta)) - ) ) -(* - | 1, 0 - | 0, 1 -> Some (Determinant.negate_phase k) - | 0, 1 -> Some (Determinant.vac 1) -*) - | _ -> None - - - -let dressing_vector ~frozen_core aux_basis f12_amplitudes ci = +let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = if Parallel.master then Printf.printf "Building matrix\n%!"; - - (* Determinants of the FCI space as a list *) - let in_dets = - DeterminantSpace.determinant_stream ci.CI.det_space - |> Util.stream_to_list + let det_space = + ci.CI.det_space in - (* Stream that generates only singly and doubly excited determinants - wrt FCI space *) - let out_dets_stream = - (* Stream that generates all determinants of FCI space *) - let s = - DeterminantSpace.fci_of_mo_basis ~frozen_core aux_basis - |> DeterminantSpace.determinant_stream - - in - (* Select only doubly excited determinants wrt FCI space *) - Stream.from (fun _ -> - try - let p12 = p12 ci.CI.det_space in - let rec result () = - let ki = Stream.next s in - match p12 ki with - | Some ki' -> Some (ki, ki') - | None -> result () - in - result () - with Stream.Failure -> None - ) - in - - let make_h_and_f alpha_list = - - let rec col_vecs_list accu_H accu_F = function - | [] -> - List.rev accu_H, - List.rev accu_F - | (ki, ki') :: rest -> - begin - let h, f = - List.map (fun kj -> - match hf_ij aux_basis kj ki with - | [ a ; b ] -> a, b - | _ -> assert false ) in_dets - |> List.split - in - let f' = - List.map (fun kj -> f_ij aux_basis kj ki') in_dets - in - let h = Vec.of_list h in - let f = Vec.of_list f in - let f' = Vec.of_list f' in - scal 0.375 f; - scal 0.125 f'; - let f = Vec.add f f' in - col_vecs_list (h::accu_H) (f::accu_F) rest - end - in - let h, f = - col_vecs_list [] [] alpha_list - in - Mat.of_col_vecs_list h, - Mat.of_col_vecs_list f - in - - let m_HF = - let batch_size = 1 + 1_000_000 / (Mat.dim1 f12_amplitudes) in - let input_stream = - Stream.from (fun i -> - let rec make_batch accu = function - | 0 -> accu - | n -> try - let alpha = Stream.next out_dets_stream in - let accu = alpha :: accu in - make_batch accu (n-1) - with Stream.Failure -> accu - in - let result = make_batch [] batch_size in - if result = [] then None else Some result - ) + + let f = + match Ds.determinants det_space with + | Ds.Arbitrary _ -> CI.create_matrix_arbitrary + | Ds.Spin _ -> CI.create_matrix_spin_computed in + f (fun deg_a deg_b ki kj -> + hf_ij_non_zero hf12_integrals deg_a deg_b ki kj + ) det_space - let iteration input = - Printf.printf ".%!"; - let m_H_aux, m_F_aux = make_h_and_f input in - let m_HF = - gemm m_H_aux m_F_aux ~transb:`T - in - gemm m_HF f12_amplitudes - in - - let result = - let x = - try [ Stream.next out_dets_stream ] - with Stream.Failure -> failwith "Auxiliary basis set does not produce any excited determinant" - in - iteration x - in - - input_stream - |> Farm.run ~ordered:false ~f:iteration - |> Stream.iter (fun hf -> - ignore @@ Mat.add result hf ~c:result ); - Printf.printf "\n"; - Parallel.broadcast (lazy result) in - if Parallel.master then Printf.printf "Done\n%!"; - Matrix.dense_of_mat m_HF - - - + Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes) let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename ?(state=1) () = - let f12 = Util.of_some @@ Simulation.f12 simulation in - let mo_num = MOBasis.size mo_basis in - -Printf.printf "Add aux basis\n%!"; - (* Add auxiliary basis set *) - let s = - let charge = Charge.to_int @@ Simulation.charge simulation - and multiplicity = Electrons.multiplicity @@ Simulation.electrons simulation - and nuclei = Simulation.nuclei simulation - in - let general_basis = - Basis.general_basis @@ Simulation.basis simulation - in - GeneralBasis.combine [ - general_basis ; GeneralBasis.read aux_basis_filename - ] - |> Basis.of_nuclei_and_general_basis nuclei - |> Simulation.make ~f12 ~charge ~multiplicity ~nuclei - in - - let aux_basis = - MOBasis.of_mo_basis s mo_basis - in - let () = - ignore @@ MOBasis.f12_ints aux_basis - in - let () = - ignore @@ MOBasis.two_e_ints aux_basis - in - let det_space = - DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num + DeterminantSpace.fci_of_mo_basis mo_basis ~frozen_core in let ci = CI.make ~n_states:state det_space in + let hf12_integrals = + HF12.make ~simulation ~mo_basis ~aux_basis_filename () + in + let ci_coef, ci_energy = let x = Lazy.force ci.eigensystem in Parallel.broadcast (lazy x) in - let e_shift = - let det = - DeterminantSpace.determinant_stream det_space - |> Stream.next - in - h_ij aux_basis det det - in - let eigensystem = lazy ( let m_H = Lazy.force ci.CI.m_H @@ -320,7 +89,7 @@ Printf.printf "Add aux basis\n%!"; let delta = (* delta_i = {% $\sum_j c_j H_{ij}$ %} *) - dressing_vector ~frozen_core aux_basis psi ci + dressing_vector ~frozen_core hf12_integrals psi ci |> Matrix.to_mat in @@ -449,14 +218,14 @@ Printf.printf "Add aux basis\n%!"; (Mat.to_col_vecs eigenvectors).(0) ) in if Parallel.master then - Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{state} +. e_shift + Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{state} +. Simulation.nuclear_repulsion simulation); if conv > threshold then iteration ~state eigenvectors else let eigenvalues = - Vec.map (fun x -> x +. e_shift) eigenvalues + Vec.map (fun x -> x +. ci.CI.e_shift) eigenvalues in eigenvectors, eigenvalues in @@ -464,6 +233,8 @@ Printf.printf "Add aux basis\n%!"; ) in - { mo_basis ; aux_basis ; det_space ; ci ; eigensystem } + { mo_basis ; det_space ; ci ; hf12_integrals ; eigensystem } + + diff --git a/MOBasis/HF12.ml b/MOBasis/HF12.ml new file mode 100644 index 0000000..10c352c --- /dev/null +++ b/MOBasis/HF12.ml @@ -0,0 +1,186 @@ +(** %{ $ \langle ij | H F | kl \rangle $ %} integrals. *) + +open Lacaml.D + +module Fis = FourIdxStorage + +type t = (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t + * (float, Bigarray.float64_elt, Bigarray.fortran_layout) Bigarray.Genarray.t + + +let make ~simulation ~mo_basis ~aux_basis_filename () = + + let f12 = Util.of_some @@ Simulation.f12 simulation in + let mo_num = MOBasis.size mo_basis in + + (* Add auxiliary basis set *) + let aux_basis = + let s = + let charge = Charge.to_int @@ Simulation.charge simulation + and multiplicity = Electrons.multiplicity @@ Simulation.electrons simulation + and nuclei = Simulation.nuclei simulation + in + let general_basis = + Basis.general_basis @@ Simulation.basis simulation + in + GeneralBasis.combine [ + general_basis ; GeneralBasis.read aux_basis_filename + ] + |> Basis.of_nuclei_and_general_basis nuclei + |> Simulation.make ~f12 ~charge ~multiplicity ~nuclei + in + MOBasis.of_mo_basis s mo_basis + in + + let aux_num = MOBasis.size aux_basis in + + (* Fire calculation of F12 and ERI *) + let f12 = + MOBasis.f12_ints aux_basis + in + let eri = + MOBasis.two_e_ints aux_basis + in + + (* Compute the integrals *) + if Parallel.master then Printf.eprintf "Computing HF12 integrals\n%!"; + + let result_s, result_o = + Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num |] , + Bigarray.Genarray.create Float64 Bigarray.fortran_layout [| mo_num ; mo_num ; mo_num ; mo_num |] + in + + + let h_s = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num aux_num aux_num in + let f_s = Bigarray.Array3.create Float64 Bigarray.fortran_layout aux_num aux_num mo_num in + let h_o = Bigarray.Array3.create Float64 Bigarray.fortran_layout mo_num aux_num aux_num in + let f_o = Bigarray.Array3.create Float64 Bigarray.fortran_layout aux_num aux_num mo_num in + let hf_s = Mat.create mo_num mo_num in + let hf_o = Mat.create mo_num mo_num in + + for a=1 to mo_num do + for b=1 to mo_num do + for i=1 to mo_num do + h_s.{i, a, b} <- 0. ; + h_o.{i, a, b} <- 0. + done + done + done; + + for k=1 to mo_num do + for b=1 to mo_num do + for a=1 to mo_num do + f_s.{a, b, k} <- 0. ; + f_o.{a, b, k} <- 0. + done + done + done; + + let task (j,l) = + + let h i a b = + h_s.{i, a, b} <- ERI.get_phys eri i j a b -. ERI.get_phys eri i j b a ; + h_o.{i, a, b} <- ERI.get_phys eri i j a b + and f a b k = + f_s.{a, b, k} <- 0.25 *. (F12.get_phys f12 a b k l -. F12.get_phys f12 a b l k) ; + f_o.{a, b, k} <- 0.375 *. F12.get_phys f12 a b k l +. 0.125 *. F12.get_phys f12 b a k l + in + + for a=mo_num+1 to aux_num do + for b=mo_num+1 to aux_num do + for i=1 to mo_num do + h i a b + done + done + done; + + for k=1 to mo_num do + for b=mo_num+1 to aux_num do + for a=mo_num+1 to aux_num do + f a b k + done + done + done; + +(* + for a=1 to mo_num do + for b=mo_num+1 to aux_num do + for i=1 to mo_num do + if i <> a then + h i a b + done + done + done; + + for k=1 to mo_num do + for b=mo_num+1 to aux_num do + for a=1 to mo_num do + if k <> a then + f a b k + done + done + done; +*) + let h_o = + Bigarray.(reshape (genarray_of_array3 h_o)) [| mo_num ; aux_num*aux_num |] + |> Bigarray.array2_of_genarray + in + let f_o = + Bigarray.(reshape (genarray_of_array3 f_o)) [| aux_num*aux_num ; mo_num |] + |> Bigarray.array2_of_genarray + in + let h_s = + Bigarray.(reshape (genarray_of_array3 h_s)) [| mo_num ; aux_num*aux_num |] + |> Bigarray.array2_of_genarray + in + let f_s = + Bigarray.(reshape (genarray_of_array3 f_s)) [| aux_num*aux_num ; mo_num |] + |> Bigarray.array2_of_genarray + in + let hf_s = gemm ~alpha:1.0 ~c:hf_s h_s f_s in + let hf_o = gemm ~alpha:1.0 ~c:hf_o h_o f_o in + hf_s, hf_o, j, l + in + + let tasks = + let rec next accu = function + | _, 0 -> accu + | 0, l -> next accu (mo_num, l-1) + | j, l -> next ((j,l) :: accu) ((j-1), l) + in + next [] (mo_num, mo_num) + |> Stream.of_list + in + + + Farm.run ~f:task ~ordered:true tasks + |> Stream.iter (fun (hf_s, hf_o, j, l) -> + (* + Printf.printf "%d %d\n" j l ; + *) + for k=1 to mo_num do + for i=1 to mo_num do + result_s.{i,k,j,l} <- hf_s.{i,k} ; + result_o.{i,k,j,l} <- hf_o.{i,k} + done + done ); + + + + (* + for l=1 to mo_num do + for k=1 to mo_num do + for j=1 to mo_num do + for i=1 to mo_num do + Printf.printf "%d %d %d %d %e\n" i j k l result.{i,j,k,l} + done + done + done + done; + Printf.printf "%!"; + *) + + Parallel.broadcast (lazy (result_s, result_o) ) + + + diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index 5fa0053..ebf2f4c 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -258,6 +258,9 @@ let get_phys t i j k l = get ~r1:{ first=i ; second=k } ~r2:{ first=j ; second=l let set_chem t i j k l value = set ~r1:{ first=i ; second=j } ~r2:{ first=k ; second=l } ~value t let set_phys t i j k l value = set ~r1:{ first=i ; second=k } ~r2:{ first=j ; second=l } ~value t +let increment_chem t i j k l value = increment ~r1:{ first=i ; second=j } ~r2:{ first=k ; second=l } ~value t +let increment_phys t i j k l value = increment ~r1:{ first=i ; second=k } ~r2:{ first=j ; second=l } ~value t + type element = (** Element for the stream *) diff --git a/Utils/FourIdxStorage.mli b/Utils/FourIdxStorage.mli index 9bddd54..aa60e49 100644 --- a/Utils/FourIdxStorage.mli +++ b/Utils/FourIdxStorage.mli @@ -38,6 +38,12 @@ val set_chem : t -> int -> int -> int -> int -> float -> unit val set_phys : t -> int -> int -> int -> int -> float -> unit (** Set an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *) +val increment_chem : t -> int -> int -> int -> int -> float -> unit +(** Increments an integral using the Chemist's convention {% $(ij|kl)$ %}. *) + +val increment_phys : t -> int -> int -> int -> int -> float -> unit +(** Increments an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *) + val get_chem_all_i : t -> j:int -> k:int -> l:int -> Vec.t (** Get all integrals in an array [a.{i} =] {% $(\cdot j|kl)$ %} . *) diff --git a/run_fci_f12.ml b/run_fci_f12.ml index d7e3390..e7104f2 100644 --- a/run_fci_f12.ml +++ b/run_fci_f12.ml @@ -87,6 +87,7 @@ let () = MOBasis.of_hartree_fock hf in + let fcif12 = F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ~state () in