From 3874dd1d951847f25d7b384e462e2e4294b467a2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 21 Mar 2019 00:44:10 +0100 Subject: [PATCH] Moved 4-idx into FourIndex.mli --- Basis/Basis.ml | 7 ++- Basis/Basis.mli | 2 + Basis/GeneralBasis.ml | 16 +++--- Basis/GeneralBasis.mli | 7 +++ CI/CI.ml | 12 ++--- CI/F12CI.ml | 102 +++++++++++++++++++++++++++++++++++++ MOBasis/MOBasis.ml | 106 ++++----------------------------------- MOBasis/MOBasis.mli | 3 ++ Utils/FourIdxStorage.ml | 100 ++++++++++++++++++++++++++++++++++++ Utils/FourIdxStorage.mli | 4 ++ run_fci_f12.ml | 33 ++++-------- 11 files changed, 257 insertions(+), 135 deletions(-) create mode 100644 CI/F12CI.ml diff --git a/Basis/Basis.ml b/Basis/Basis.ml index dec6c92..cae8fd8 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -1,8 +1,9 @@ type t = { - size : int; + size : int; contracted_shells : ContractedShell.t array ; atomic_shells : AtomicShell.t array lazy_t; + general_basis : GeneralBasis.t ; } module As = AtomicShell @@ -10,6 +11,7 @@ module Cs = ContractedShell module Gb = GeneralBasis module Ps = PrimitiveShell +let general_basis t = t.general_basis (** Returns an array of the basis set per atom *) let of_nuclei_and_general_basis nucl bas = @@ -49,7 +51,8 @@ let of_nuclei_and_general_basis nucl bas = |> List.sort (fun x y -> compare (As.index x) (As.index y)) |> Array.of_list ) in - { contracted_shells ; atomic_shells ; size = !index_ } + { contracted_shells ; atomic_shells ; size = !index_; + general_basis = bas } let size x = x.size diff --git a/Basis/Basis.mli b/Basis/Basis.mli index b5e43d2..9402864 100644 --- a/Basis/Basis.mli +++ b/Basis/Basis.mli @@ -33,6 +33,8 @@ val atomic_shells : t -> AtomicShell.t array val contracted_shells : t -> ContractedShell.t array (** Returns all the contracted basis functions. *) +val general_basis : t -> GeneralBasis.t +(** Returns the [!GeneralBasis] that was used to build the current basis. *) (** {2 Printers} *) diff --git a/Basis/GeneralBasis.ml b/Basis/GeneralBasis.ml index bf2c5d5..c5e186f 100644 --- a/Basis/GeneralBasis.ml +++ b/Basis/GeneralBasis.ml @@ -78,18 +78,20 @@ let read filename = loop () -let read_many filenames = +let combine basis_list = let h = Hashtbl.create 63 in - let l = - List.map read filenames - |> List.concat - in - List.iter (fun (k,v) -> + List.concat basis_list + |> List.iter (fun (k,v) -> let l = Hashtbl.find_all h k in Hashtbl.replace h k (Array.concat (v::l) ) - ) l; + ) ; Hashtbl.fold (fun k v accu -> (k, v)::accu) h [] + +let read_many filenames = + List.map read filenames + |> combine + let string_of_primitive ?id prim = match id with diff --git a/Basis/GeneralBasis.mli b/Basis/GeneralBasis.mli index 45421c6..138022a 100644 --- a/Basis/GeneralBasis.mli +++ b/Basis/GeneralBasis.mli @@ -51,6 +51,12 @@ val read : string -> t *) +val combine : t list -> t +(** Combines the basis functions of each element among different basis sets. + Used to complement a basis with an auxiliary basis set. +*) + + val read_many : string list -> t (** Reads multiple basis set files and return an association list where the key is an {!Element.t} and the value is the parsed basis set. @@ -58,6 +64,7 @@ val read_many : string list -> t *) + val read_element : in_channel -> element_basis option (** Reads an element from the input channel [ic]. For example, {[ diff --git a/CI/CI.ml b/CI/CI.ml index cad76d5..9133690 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -592,12 +592,9 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } same_spin +. opposite_spin in - 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) + Util.stream_range 0 (Array.length psi0 - 1) + |> Farm.run ~ordered:true ~f:det_contribution + |> Util.stream_to_list let pt2_en ci = @@ -662,6 +659,7 @@ let pt2_en ci = second_order_sum ci list_holes list_particles i_o1_alfa i_o1_alfa w_alfa + |> List.fold_left (+.) 0. @@ -692,6 +690,7 @@ let pt2_mp ci = second_order_sum ci list_holes list_particles i_o1_alfa i_o1_alfa w_alfa + |> List.fold_left (+.) 0. let variance ci = @@ -711,5 +710,6 @@ let variance ci = second_order_sum ci list_holes list_particles i_o1_alfa i_o1_alfa w_alfa + |> List.fold_left (+.) 0. diff --git a/CI/F12CI.ml b/CI/F12CI.ml new file mode 100644 index 0000000..d02a831 --- /dev/null +++ b/CI/F12CI.ml @@ -0,0 +1,102 @@ +open Lacaml.D + +type t = + { + mo_basis : MOBasis.t ; + aux_basis : MOBasis.t ; + det_space : DeterminantSpace.t ; + ci : CI.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 f12_integrals mo_basis = + let two_e_ints = MOBasis.two_e_ints mo_basis in + ( (fun i j _ -> 0.), + (fun i j k l s s' -> + if s' = Spin.other s then + ERI.get_phys two_e_ints i j k l + else + (ERI.get_phys two_e_ints i j k l) -. + (ERI.get_phys two_e_ints i j l k) + ) ) + + + +let h_ij mo_basis ki kj = + let integrals = + List.map (fun f -> f mo_basis) + [ CI.h_integrals ] + 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 + |> List.hd + + +let dressing_vector ci = + + let mo_basis = DeterminantSpace.mo_basis ci.CI.det_space in + + let i_o1_alfa = h_ij mo_basis in + let alfa_o2_i = f_ij mo_basis in + + let w_alfa _ _ = 1. in + + let mo_class = CI.mo_class ci in + let list_holes = List.concat + [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ] + and list_particles = MOClass.auxiliary_mos mo_class + in + + CI.second_order_sum ci list_holes list_particles + i_o1_alfa alfa_o2_i w_alfa + |> Vec.of_list + + + + + + +let make ~simulation ?(frozen_core=true) ~mo_basis ~aux_basis_filename () = + + let mo_num = MOBasis.size mo_basis in + + (* 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 ~charge ~multiplicity ~nuclei + in + + let aux_basis = + MOBasis.of_mo_basis s mo_basis + in + + let det_space = + DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num + in + + let ci = CI.make det_space in + { mo_basis ; aux_basis ; det_space ; ci } + + diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml index 4e2e8e7..796bbbe 100644 --- a/MOBasis/MOBasis.ml +++ b/MOBasis/MOBasis.ml @@ -20,6 +20,7 @@ type t = mo_coef : Mat.t; (* Matrix of the MO coefficients in the AO basis *) eN_ints : NucInt.t lazy_t; (* Electron-nucleus potential integrals *) ee_ints : ERI.t lazy_t; (* Electron-electron potential integrals *) + f12_ints : F12.t lazy_t; (* F12 integrals *) kin_ints : KinInt.t lazy_t; (* Kinetic energy integrals *) one_e_ints : Mat.t lazy_t; (* Kinetic energy integrals *) } @@ -37,6 +38,7 @@ let eN_ints t = Lazy.force t.eN_ints let ee_ints t = Lazy.force t.ee_ints let kin_ints t = Lazy.force t.kin_ints let two_e_ints t = Lazy.force t.ee_ints +let f12_ints t = Lazy.force t.f12_ints let one_e_ints t = Lazy.force t.one_e_ints let mo_energies t = @@ -64,99 +66,6 @@ let ao_matrix_of_mo_matrix ~mo_coef ~ao_overlap mo_matrix = x_o_xt ~x:sc ~o:mo_matrix -let four_index_transform ~mo_coef eri_ao = - - let ao_num = Mat.dim1 mo_coef in - let mo_num = Mat.dim2 mo_coef in - let eri_mo = - ERI.create ~size:mo_num `Dense - in - - let mo_num_2 = mo_num * mo_num in - let ao_num_2 = ao_num * ao_num in - let ao_mo_num = ao_num * mo_num in - - let range_mo = list_range 1 mo_num in - let range_ao = list_range 1 ao_num in - - let u = Mat.create mo_num_2 mo_num - and o = Mat.create ao_num ao_num_2 - and p = Mat.create ao_num_2 mo_num - and q = Mat.create ao_mo_num mo_num - in - - if Parallel.master then Printf.eprintf "4-idx transformation \n%!"; - - let task delta = - Mat.fill u 0.; - - List.iter (fun l -> - if abs_float mo_coef.{l,delta} > epsilon then - begin - let jk = ref 0 in - List.iter (fun k -> - List.iter (fun j -> - incr jk; - ERI.get_chem_all_i eri_ao ~j ~k ~l - |> Array.iteri (fun i x -> o.{i+1,!jk} <- x) - ) range_ao - ) range_ao; - (* o_i_jk *) - - let p = - gemm ~transa:`T ~c:p o mo_coef - (* p_jk_alpha = \sum_i o_i_jk c_i_alpha *) - in - let p' = - Bigarray.reshape_2 (Bigarray.genarray_of_array2 p) ao_num ao_mo_num - (* p_j_kalpha *) - in - - let q = - gemm ~transa:`T ~c:q p' mo_coef - (* q_kalpha_beta = \sum_j p_j_kalpha c_j_beta *) - in - let q' = - Bigarray.reshape_2 (Bigarray.genarray_of_array2 q) ao_num mo_num_2 - (* q_k_alphabeta = \sum_j p_j_kalpha c_j_beta *) - in - - ignore @@ - gemm ~transa:`T ~beta:1. ~alpha:mo_coef.{l,delta} ~c:u q' mo_coef ; - (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *) - end - ) range_ao; - let u = - Bigarray.reshape - (Bigarray.genarray_of_array2 u) - [| mo_num ; mo_num ; mo_num |] - |> Bigarray.array3_of_genarray - in - let result = ref [] in - List.iter (fun gamma -> - List.iter (fun beta -> - List.iter (fun alpha -> - let x = u.{alpha,beta,gamma} in - if x <> 0. then - result := (alpha, beta, gamma, delta, x) :: !result; - ) (list_range 1 beta) - ) range_mo - ) (list_range 1 delta); - Array.of_list !result - in - - let n = ref 0 in - Stream.of_list range_mo - |> Farm.run ~f:task ~ordered:false - |> Stream.iter (fun l -> - if Parallel.master then (Printf.eprintf "\r%d / %d%!" !n mo_num; incr n); - Array.iter (fun (alpha, beta, gamma, delta, x) -> - ERI.set_chem eri_mo alpha beta gamma delta x) l); - - if Parallel.master then Printf.eprintf "\n"; - Parallel.broadcast (lazy eri_mo) - - let make ~simulation ~mo_type ~mo_occupation ~mo_coef () = let ao_basis = Si.ao_basis simulation @@ -164,7 +73,7 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef () = let eN_ints = lazy ( AOBasis.eN_ints ao_basis |> NucInt.matrix - |> mo_matrix_of_ao_matrix ~mo_coef + |> mo_matrix_of_ao_matrix ~mo_coef |> NucInt.of_matrix ) and kin_ints = lazy ( @@ -175,7 +84,11 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef () = ) and ee_ints = lazy ( AOBasis.ee_ints ao_basis - |> four_index_transform ~mo_coef + |> ERI.four_index_transform mo_coef + ) + and f12_ints = lazy ( + AOBasis.f12_ints ao_basis + |> F12.four_index_transform mo_coef ) in let one_e_ints = lazy ( @@ -183,7 +96,8 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef () = (KinInt.matrix @@ Lazy.force kin_ints) ) in { simulation ; mo_type ; mo_occupation ; mo_coef ; - eN_ints ; ee_ints ; kin_ints ; one_e_ints } + eN_ints ; ee_ints ; kin_ints ; one_e_ints ; + f12_ints } diff --git a/MOBasis/MOBasis.mli b/MOBasis/MOBasis.mli index 25ef20c..0c6e055 100644 --- a/MOBasis/MOBasis.mli +++ b/MOBasis/MOBasis.mli @@ -46,6 +46,9 @@ val one_e_ints : t -> Mat.t val two_e_ints : t -> ERI.t (** Electron-electron repulsion integrals *) +val f12_ints : t -> F12.t +(** F12 integrals *) + val size : t -> int (** Number of molecular orbitals in the basis *) diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index d9a89aa..376d511 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -1,4 +1,6 @@ open Util +open Lacaml.D +open Constants let max_index = 1 lsl 14 @@ -260,3 +262,101 @@ let to_list data = in append [] + + + + +let four_index_transform coef source = + + let ao_num = Mat.dim1 coef in + let mo_num = Mat.dim2 coef in + let destination = + match source.four_index with + | Dense _ -> create ~size:mo_num `Dense + | Sparse _ -> create ~size:mo_num `Sparse + in + + let mo_num_2 = mo_num * mo_num in + let ao_num_2 = ao_num * ao_num in + let ao_mo_num = ao_num * mo_num in + + let range_mo = list_range 1 mo_num in + let range_ao = list_range 1 ao_num in + + let u = Mat.create mo_num_2 mo_num + and o = Mat.create ao_num ao_num_2 + and p = Mat.create ao_num_2 mo_num + and q = Mat.create ao_mo_num mo_num + in + + if Parallel.master then Printf.eprintf "4-idx transformation \n%!"; + + let task delta = + Mat.fill u 0.; + + List.iter (fun l -> + if abs_float coef.{l,delta} > epsilon then + begin + let jk = ref 0 in + List.iter (fun k -> + List.iter (fun j -> + incr jk; + get_chem_all_i source ~j ~k ~l + |> Array.iteri (fun i x -> o.{i+1,!jk} <- x) + ) range_ao + ) range_ao; + (* o_i_jk *) + + let p = + gemm ~transa:`T ~c:p o coef + (* p_jk_alpha = \sum_i o_i_jk c_i_alpha *) + in + let p' = + Bigarray.reshape_2 (Bigarray.genarray_of_array2 p) ao_num ao_mo_num + (* p_j_kalpha *) + in + + let q = + gemm ~transa:`T ~c:q p' coef + (* q_kalpha_beta = \sum_j p_j_kalpha c_j_beta *) + in + let q' = + Bigarray.reshape_2 (Bigarray.genarray_of_array2 q) ao_num mo_num_2 + (* q_k_alphabeta = \sum_j p_j_kalpha c_j_beta *) + in + + ignore @@ + gemm ~transa:`T ~beta:1. ~alpha:coef.{l,delta} ~c:u q' coef ; + (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *) + end + ) range_ao; + let u = + Bigarray.reshape + (Bigarray.genarray_of_array2 u) + [| mo_num ; mo_num ; mo_num |] + |> Bigarray.array3_of_genarray + in + let result = ref [] in + List.iter (fun gamma -> + List.iter (fun beta -> + List.iter (fun alpha -> + let x = u.{alpha,beta,gamma} in + if x <> 0. then + result := (alpha, beta, gamma, delta, x) :: !result; + ) (list_range 1 beta) + ) range_mo + ) (list_range 1 delta); + Array.of_list !result + in + + let n = ref 0 in + Stream.of_list range_mo + |> Farm.run ~f:task ~ordered:false + |> Stream.iter (fun l -> + if Parallel.master then (Printf.eprintf "\r%d / %d%!" !n mo_num; incr n); + Array.iter (fun (alpha, beta, gamma, delta, x) -> + set_chem destination alpha beta gamma delta x) l); + + if Parallel.master then Printf.eprintf "\n"; + Parallel.broadcast (lazy destination) + diff --git a/Utils/FourIdxStorage.mli b/Utils/FourIdxStorage.mli index 3129200..28beef6 100644 --- a/Utils/FourIdxStorage.mli +++ b/Utils/FourIdxStorage.mli @@ -7,6 +7,8 @@ There are two kinds of ordering of indices: *) +open Lacaml.D + type t type element = (** Element for the stream *) @@ -54,6 +56,8 @@ val to_stream : t -> element Stream.t val to_list : t -> element list (** Retrun the data structure as a list. *) +val four_index_transform : Mat.t -> t -> t +(** Perform a four-index transformation *) (** {2 I/O} *) diff --git a/run_fci_f12.ml b/run_fci_f12.ml index eae08d3..ddc20f1 100644 --- a/run_fci_f12.ml +++ b/run_fci_f12.ml @@ -2,7 +2,7 @@ let () = let open Command_line in begin set_header_doc (Sys.argv.(0) ^ " - QCaml command"); - set_description_doc "Runs a Hartree-Fock calculation"; + set_description_doc "Runs a F12-Full CI calculation"; set_specs [ { short='b' ; long="basis" ; opt=Mandatory; arg=With_arg ""; @@ -29,7 +29,7 @@ let () = (* Handle options *) let basis_file = Util.of_some @@ Command_line.get "basis" in - let aux_basis_file = Util.of_some @@ Command_line.get "aux_basis" in + let aux_basis_filename = Util.of_some @@ Command_line.get "aux_basis" in let nuclei_file = Util.of_some @@ Command_line.get "xyz" in @@ -50,37 +50,22 @@ let () = else Printing.ppf_dev_null in - let s' = + let simulation = Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file in - let hf = HartreeFock.make s' in + let hf = HartreeFock.make simulation in Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf; - let mos' = + let mo_basis = MOBasis.of_hartree_fock hf in - let mo_num = - MOBasis.size mos' + let fcif12 = + F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename () in - - - let s = - Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file - ~aux_basis_filenames:[aux_basis_file] basis_file - in - - let mos = - MOBasis.of_mo_basis s mos' - in - - let space = - DeterminantSpace.fci_f12_of_mo_basis mos ~frozen_core:false mo_num - in - - let ci = CI.make space in - Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); + let ci = F12CI.ci fcif12 in + Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion simulation); (* let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in Util.list_range 1 (DeterminantSpace.size space)