diff --git a/Basis/AOBasis.ml b/Basis/AOBasis.ml index 3779f91..c4700cb 100644 --- a/Basis/AOBasis.ml +++ b/Basis/AOBasis.ml @@ -47,14 +47,8 @@ let make ~cartesian ~basis ?f12 nuclei = ) in let f12_ints = lazy ( - let f12 = - match f12 with - | Some f12 -> f12 - | None -> failwith "Missing f12 factor" - in - F12.of_basis f12 basis - ) - in + F12.of_basis basis + ) in { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; f12_ints ; cartesian ; diff --git a/Basis/AOBasis.mli b/Basis/AOBasis.mli index f92d421..0ad81c5 100644 --- a/Basis/AOBasis.mli +++ b/Basis/AOBasis.mli @@ -31,7 +31,7 @@ val cartesian : t -> bool (** {1 Creators} *) -val make : cartesian:bool -> basis:Basis.t -> ?f12:F12.f12_factor -> Nuclei.t -> t +val make : cartesian:bool -> basis:Basis.t -> ?f12:F12factor.t -> Nuclei.t -> t (** Creates the data structure for atomic orbitals from a {Basis.t} and the molecular geometry {Nuclei.t} *) diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 9d742ed..1462d74 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -3,9 +3,10 @@ open Constants open Util +module Csp = ContractedShellPair +module Cspc = ContractedShellPairCouple - -module Zm = struct +module T = struct let name = "Electron repulsion integrals" @@ -35,8 +36,20 @@ module Zm = struct aux f 0 maxm; result + let class_of_contracted_shell_pair_couple shell_pair_couple = + let shell_p = Cspc.shell_pair_p shell_pair_couple + and shell_q = Cspc.shell_pair_q shell_pair_couple + in + if Array.length (Csp.shell_pairs shell_p) + + (Array.length (Csp.shell_pairs shell_q)) < 4 then + TwoElectronRR.contracted_class_shell_pair_couple + ~zero_m shell_pair_couple + else + TwoElectronRRVectorized.contracted_class_shell_pairs + ~zero_m shell_p shell_q + end -module M = TwoElectronIntegralsNonSeparable.Make(Zm) +module M = TwoElectronIntegrals.Make(T) include M diff --git a/Basis/F12.ml b/Basis/F12.ml index 3f0c79f..c24b64b 100644 --- a/Basis/F12.ml +++ b/Basis/F12.ml @@ -5,397 +5,23 @@ open Constants let cutoff = integrals_cutoff -module Bs = Basis -module Cs = ContractedShell -module Csp = ContractedShellPair -module Cspc = ContractedShellPairCouple -module Fis = FourIdxStorage +module T = struct - -include FourIdxStorage - -type f12_factor = -{ - expo_s : float ; - coef_g : float array; - expo_sg : float array; - expo_sg_inv : float array; -} - -let make_gaussian_corr_factor expo_s coef_g expo_sg = - let expo_sg_inv = - Array.map (fun x -> 1. /. (x *. expo_s *. expo_s)) expo_sg - in - { - expo_s ; coef_g ; expo_sg ; expo_sg_inv - } - - -(* -1/expo_s *. exp (-expo_s r) *) -let gaussian_geminal expo_s = - let coef_g = - [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |] - |> Array.map (fun x -> -. x /. expo_s) - and expo_sg = - [| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |] - in - make_gaussian_corr_factor expo_s coef_g expo_sg - - - -(* exp (-expo_s r) *) -let simple_gaussian_geminal expo_s = - let coef_g = - [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |] - and expo_sg = - [| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |] - in - make_gaussian_corr_factor expo_s coef_g expo_sg - - - -(** r12 * exp ( -expo_s * r) *) -let gaussian_geminal_times_r12 expo_s = - let coef_g = - [| 0.2454 ; 0.2938 ; 0.1815 ; 0.11281 ; 0.07502 ; 0.05280 |] - and expo_sg = - [| 0.1824 ; 0.7118; 2.252 ; 6.474 ; 19.66 ; 77.92 |] - in make_gaussian_corr_factor expo_s coef_g expo_sg - - -(* exp (-expo_s r) *) -let simple_gaussian_geminal' expo_s = - let coef_g = - [| - -3.4793465193721626604883567779324948787689208984375 ; - -0.00571703486454788484955047422886309504974633455276489257812 ; - 4.14878218728681513738365538301877677440643310546875 ; - 0.202874298181392742623785352407139725983142852783203125 ; - 0.0819187742387294803858566183407674543559551239013671875 ; - 0.04225945671351955673644695821167260874062776565551757812 ; - |] - and expo_sg = - [| - 0.63172472556807146570889699432882480323314666748046875; - 26.3759196683467962429858744144439697265625; - 0.63172102793029016876147352377302013337612152099609375; - 7.08429025944207335641067402320913970470428466796875; - 42.4442841447001910637482069432735443115234375; - 391.44036073596890901171718724071979522705078125 ; - |] - in make_gaussian_corr_factor expo_s coef_g expo_sg - - - - -let one_over_r = - - let coef_g = [| - 841.88478132 ; 70.590185207 ; 18.3616020768 ; 7.2608642093 ; - 3.57483416444 ; 2.01376031082 ; 1.24216542801 ; 0.81754348620 ; - 0.564546514023 ; 0.404228610699 ; 0.297458536575 ; 0.223321219537 ; - 0.169933732064 ; 0.130190978230 ; 0.099652303426 ; 0.075428246546 ; - 0.0555635614051 ; 0.0386791283055 ; 0.0237550435652 ; 0.010006278387 ; - |] - and expo_sg_inv = - [| 84135.654509 ; 2971.58727634 ; 474.716025959 ; 130.676724560 ; - 47.3938388887 ; 20.2078651631 ; 9.5411021938 ; 4.8109546955 ; - 2.52795733067 ; 1.35894103210 ; 0.73586710268 ; 0.39557629706 ; - 0.20785895177 ; 0.104809693858 ; 0.049485682527 ; 0.021099788990 ; - 0.007652472186 ; 0.0021065225215 ; 0.0003365204879 ; 0.0000118855674 |] - in make_gaussian_corr_factor 1.0 coef_g expo_sg_inv - - - - - - -module Zero_m = struct let name = "F12" -end -let class_of_contracted_shell_pair_couple f12 shell_pair_couple = + let f12_factor = F12factor.gaussian_geminal 1.0 + + let class_of_contracted_shell_pair_couple shell_pair_couple = + let g = f12_factor.F12factor.gaussian in F12RR.contracted_class_shell_pair_couple - f12.expo_sg_inv f12.coef_g shell_pair_couple - - -let filter_contracted_shell_pairs f12 ?(cutoff=integrals_cutoff) shell_pairs = - List.map (fun pair -> - match Cspc.make ~cutoff pair pair with - | Some cspc -> - let cls = class_of_contracted_shell_pair_couple f12 cspc in - (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) - (* TODO \sum_k |coef_k * integral_k| *) - | None -> (pair, -1.) - ) shell_pairs - |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) - |> List.map fst - - -(* TODO - let filter_contracted_shell_pair_couples - ?(cutoff=integrals_cutoff) shell_pair_couples = - List.map (fun pair -> - let cls = - class_of_contracted_shell_pairs pair pair - in - (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) - ) shell_pairs - |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) - |> List.map fst - -*) - - -let store_class basis ?(cutoff=integrals_cutoff) data f12 contracted_shell_pair_couple cls = - let to_powers x = - let open Zkey in - match to_powers x with - | Three x -> x - | _ -> assert false - in - - let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple - and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple - in - - (* - let lambda_inv = -. 1. /. f12.expo_s in - *) - Array.iteri (fun i_c powers_i -> - let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in - let xi = to_powers powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in - let xj = to_powers powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in - let xk = to_powers powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in - let xl = to_powers powers_l in - let key = Zkey.of_powers_twelve xi xj xk xl in - let value = Zmap.find cls key in - (* - lambda_inv *. (s.{i_c,j_c} *. s.{k_c,l_c} -. value) - lambda_inv *. value - *) - value - |> set_chem data i_c j_c k_c l_c - ) (Cs.zkey_array (Csp.shell_b shell_q)) - ) (Cs.zkey_array (Csp.shell_a shell_q)) - ) (Cs.zkey_array (Csp.shell_b shell_p)) - ) (Cs.zkey_array (Csp.shell_a shell_p)) - - - - - - -let of_basis_serial f12 basis = - - let n = Bs.size basis - and shell = Bs.contracted_shells basis - in - - let eri_array = - Fis.create ~size:n `Dense -(* - Fis.create ~size:n `Sparse -*) - in - - let t0 = Unix.gettimeofday () in - - let shell_pairs = - Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs f12 ~cutoff - in - - Printf.printf "%d significant shell pairs computed in %f seconds\n" - (List.length shell_pairs) (Unix.gettimeofday () -. t0); - - - let t0 = Unix.gettimeofday () in - let ishell = ref 0 in - - List.iter (fun shell_p -> - let () = - if (Cs.index (Csp.shell_a shell_p) > !ishell) then - (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) - in - - let sp = - Csp.shell_pairs shell_p - in - - try - List.iter (fun shell_q -> - let () = - if Cs.index (Csp.shell_a shell_q) > - Cs.index (Csp.shell_a shell_p) then - raise Exit - in - let sq = Csp.shell_pairs shell_q in - let cspc = - if Array.length sp < Array.length sq then - Cspc.make ~cutoff shell_p shell_q - else - Cspc.make ~cutoff shell_q shell_p - in - - match cspc with - | Some cspc -> - let cls = - class_of_contracted_shell_pair_couple f12 cspc - in - store_class basis ~cutoff eri_array f12 cspc cls - | None -> () - ) shell_pairs - with Exit -> () - ) shell_pairs ; - Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0); - eri_array - - - - - - - -(* Parallel functions *) - - - -let of_basis_parallel f12 basis = - - let n = Bs.size basis - and shell = Bs.contracted_shells basis - in - - let store_class_parallel - ?(cutoff=integrals_cutoff) contracted_shell_pair_couple cls = - let to_powers x = - let open Zkey in - match to_powers x with - | Three x -> x - | _ -> assert false - in - - let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple - and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple - in - - let result = ref [] in - Array.iteri (fun i_c powers_i -> - let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in - let xi = to_powers powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in - let xj = to_powers powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in - let xk = to_powers powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in - let xl = to_powers powers_l in - let key = Zkey.of_powers_twelve xi xj xk xl in - let value = Zmap.find cls key in - result := (i_c, j_c, k_c, l_c, value) :: !result - ) (Cs.zkey_array (Csp.shell_b shell_q)) - ) (Cs.zkey_array (Csp.shell_a shell_q)) - ) (Cs.zkey_array (Csp.shell_b shell_p)) - ) (Cs.zkey_array (Csp.shell_a shell_p)); - !result - in - - - - let t0 = Unix.gettimeofday () in - - let shell_pairs = - Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs f12 ~cutoff - in - - if Parallel.master then - Printf.printf "%d significant shell pairs computed in %f seconds\n" - (List.length shell_pairs) (Unix.gettimeofday () -. t0); - - - let t0 = Unix.gettimeofday () in - let ishell = ref max_int in - - let input_stream = Stream.of_list (List.rev shell_pairs) in - - let f shell_p = - let () = - if Parallel.rank < 2 && Cs.index (Csp.shell_a shell_p) < !ishell then - (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) - in - - let sp = - Csp.shell_pairs shell_p - in - - let result = ref [] in - try - List.iter (fun shell_q -> - let () = - if Cs.index (Csp.shell_a shell_q) > - Cs.index (Csp.shell_a shell_p) then - raise Exit - in - let sq = Csp.shell_pairs shell_q in - let cspc = - if Array.length sp < Array.length sq then - Cspc.make ~cutoff shell_p shell_q - else - Cspc.make ~cutoff shell_q shell_p - in - - match cspc with - | Some cspc -> - let cls = - class_of_contracted_shell_pair_couple f12 cspc - in - result := (store_class_parallel ~cutoff cspc cls) :: !result; - | None -> () - ) shell_pairs; - raise Exit - with Exit -> List.concat !result |> Array.of_list - in - - let eri_array = - if Parallel.master then - Fis.create ~size:n `Dense - else - Fis.create ~size:n `Dense - in - - Farm.run ~ordered:false ~f input_stream - |> Stream.iter (fun l -> - Array.iter (fun (i_c,j_c,k_c,l_c,value) -> - (* - lambda_inv *. (s.{i_c,j_c} *. s.{k_c,l_c} -. value) - lambda_inv *. value - *) - value - |> set_chem eri_array i_c j_c k_c l_c) l); - - if Parallel.master then - Printf.printf - "Computed %s Integrals in parallel in %f seconds\n%!" Zero_m.name (Unix.gettimeofday () -. t0); - Fis.broadcast eri_array - - - -let of_basis = - match Parallel.size with - | 1 -> of_basis_serial - | _ -> of_basis_parallel + g.GaussianOperator.expo_sg_inv + g.GaussianOperator.coef_g + shell_pair_couple +end +module M = TwoElectronIntegrals.Make(T) +include M diff --git a/Basis/F12factor.ml b/Basis/F12factor.ml new file mode 100644 index 0000000..0156399 --- /dev/null +++ b/Basis/F12factor.ml @@ -0,0 +1,78 @@ +(** Type for f12 correlation factors *) + +open Constants + +type t = + { + expo_s : float ; + gaussian : GaussianOperator.t + } + + +let make_gaussian_corr_factor expo_s coef_g expo_sg = + let expo_sg = + Array.map (fun x -> x *. expo_s *. expo_s) expo_sg + in + let gaussian = GaussianOperator.make coef_g expo_sg in + { expo_s ; gaussian } + + +(* -1/expo_s *. exp (-expo_s r) *) +let gaussian_geminal expo_s = + let coef_g = + [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |] + |> Array.map (fun x -> -. x /. expo_s) + and expo_sg = + [| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |] + in + make_gaussian_corr_factor expo_s coef_g expo_sg + + + +(* exp (-expo_s r) *) +let simple_gaussian_geminal expo_s = + let coef_g = + [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |] + and expo_sg = + [| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |] + in + make_gaussian_corr_factor expo_s coef_g expo_sg + + + +(** r12 * exp ( -expo_s * r) *) +let gaussian_geminal_times_r12 expo_s = + let coef_g = + [| 0.2454 ; 0.2938 ; 0.1815 ; 0.11281 ; 0.07502 ; 0.05280 |] + and expo_sg = + [| 0.1824 ; 0.7118; 2.252 ; 6.474 ; 19.66 ; 77.92 |] + in make_gaussian_corr_factor expo_s coef_g expo_sg + + +(* exp (-expo_s r) *) +let simple_gaussian_geminal' expo_s = + let coef_g = + [| + -3.4793465193721626604883567779324948787689208984375 ; + -0.00571703486454788484955047422886309504974633455276489257812 ; + 4.14878218728681513738365538301877677440643310546875 ; + 0.202874298181392742623785352407139725983142852783203125 ; + 0.0819187742387294803858566183407674543559551239013671875 ; + 0.04225945671351955673644695821167260874062776565551757812 ; + |] + and expo_sg = + [| + 0.63172472556807146570889699432882480323314666748046875; + 26.3759196683467962429858744144439697265625; + 0.63172102793029016876147352377302013337612152099609375; + 7.08429025944207335641067402320913970470428466796875; + 42.4442841447001910637482069432735443115234375; + 391.44036073596890901171718724071979522705078125 ; + |] + in make_gaussian_corr_factor expo_s coef_g expo_sg + + + + + + diff --git a/Basis/GaussianOperator.ml b/Basis/GaussianOperator.ml new file mode 100644 index 0000000..01a1242 --- /dev/null +++ b/Basis/GaussianOperator.ml @@ -0,0 +1,39 @@ +(** Representation for two-electron operators expressed in a Gaussian basis set. *) + +open Constants + +type t = +{ + coef_g : float array; + expo_sg : float array; + expo_sg_inv : float array; +} + +let make coef_g expo_sg = + let expo_sg_inv = + Array.map (fun x -> 1. /. x ) expo_sg + in + { coef_g ; expo_sg ; expo_sg_inv } + + +let one_over_r = + + let coef_g = [| + 841.88478132 ; 70.590185207 ; 18.3616020768 ; 7.2608642093 ; + 3.57483416444 ; 2.01376031082 ; 1.24216542801 ; 0.81754348620 ; + 0.564546514023 ; 0.404228610699 ; 0.297458536575 ; 0.223321219537 ; + 0.169933732064 ; 0.130190978230 ; 0.099652303426 ; 0.075428246546 ; + 0.0555635614051 ; 0.0386791283055 ; 0.0237550435652 ; 0.010006278387 ; + |] + and expo_sg = + [| 84135.654509 ; 2971.58727634 ; 474.716025959 ; 130.676724560 ; + 47.3938388887 ; 20.2078651631 ; 9.5411021938 ; 4.8109546955 ; + 2.52795733067 ; 1.35894103210 ; 0.73586710268 ; 0.39557629706 ; + 0.20785895177 ; 0.104809693858 ; 0.049485682527 ; 0.021099788990 ; + 0.007652472186 ; 0.0021065225215 ; 0.0003365204879 ; 0.0000118855674 |] + in make coef_g expo_sg + + + + + diff --git a/Basis/TwoElectronIntegralsNonSeparable.ml b/Basis/TwoElectronIntegrals.ml similarity index 90% rename from Basis/TwoElectronIntegralsNonSeparable.ml rename to Basis/TwoElectronIntegrals.ml index 83995ae..7d4530d 100644 --- a/Basis/TwoElectronIntegralsNonSeparable.ml +++ b/Basis/TwoElectronIntegrals.ml @@ -1,5 +1,4 @@ -(** Two electron integral functor for operators that are not separable among %{ $(x,y,z)$ %}. - It is parameterized by the [zero_m] function. +(** Two electron integrals *) open Constants @@ -11,30 +10,18 @@ module Csp = ContractedShellPair module Cspc = ContractedShellPairCouple module Fis = FourIdxStorage -module type Zero_mType = +module type TwoEI_structure = sig val name : string - val zero_m : Zero_m_parameters.t -> float array + val class_of_contracted_shell_pair_couple : ContractedShellPairCouple.t -> float Zmap.t end -module Make(Zero_m : Zero_mType) = struct +module Make(T : TwoEI_structure) = struct include FourIdxStorage - let zero_m = Zero_m.zero_m - - let class_of_contracted_shell_pair_couple shell_pair_couple = - let shell_p = Cspc.shell_pair_p shell_pair_couple - and shell_q = Cspc.shell_pair_q shell_pair_couple - in - if Array.length (Csp.shell_pairs shell_p) + - (Array.length (Csp.shell_pairs shell_q)) < 4 then - TwoElectronRR.contracted_class_shell_pair_couple - ~zero_m shell_pair_couple - else - TwoElectronRRVectorized.contracted_class_shell_pairs - ~zero_m shell_p shell_q + let class_of_contracted_shell_pair_couple = T.class_of_contracted_shell_pair_couple let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs = List.map (fun pair -> @@ -283,6 +270,7 @@ module Make(Zero_m : Zero_mType) = struct else Fis.create ~size:n `Dense in + Farm.run ~ordered:true ~f input_stream |> Stream.iter (fun l -> Array.iter (fun (i_c,j_c,k_c,l_c,value) -> @@ -290,7 +278,7 @@ module Make(Zero_m : Zero_mType) = struct if Parallel.master then Printf.printf - "Computed %s Integrals in parallel in %f seconds\n%!" Zero_m.name (Unix.gettimeofday () -. t0); + "Computed %s Integrals in parallel in %f seconds\n%!" T.name (Unix.gettimeofday () -. t0); Fis.broadcast eri_array diff --git a/Basis/TwoElectronIntegralsNonSeparable.mli b/Basis/TwoElectronIntegrals.mli similarity index 58% rename from Basis/TwoElectronIntegralsNonSeparable.mli rename to Basis/TwoElectronIntegrals.mli index 43ecb24..1c8e957 100644 --- a/Basis/TwoElectronIntegralsNonSeparable.mli +++ b/Basis/TwoElectronIntegrals.mli @@ -6,27 +6,20 @@ *) -module type Zero_mType = +module type TwoEI_structure = sig val name : string (** Name of the kind of integrals, for printing purposes. *) -val zero_m : Zero_m_parameters.t -> float array -(** The returned float array contains all the {% $(00|00)^m$ %} values, where - [m] is the index of the array. - - - [maxm] : Maximum total angular momentum - - [expo_pq_inv] : {% $1/p + 1/q$ %} where {% $p$ %} and {% $q$ %} are the - exponents of {% $\phi_p$ %} and {% $\phi_q$ %} - - [norm_pq_sq] : square of the distance between the centers of - {% $\phi_p$ %} and {% $\phi_q$ %} - +val class_of_contracted_shell_pair_couple : ContractedShellPairCouple.t -> float Zmap.t +(** Returns an integral class from a couple of contracted shells. + The results is stored in a Zmap. *) - end -module Make : functor (Zero_m : Zero_mType) -> + +module Make : functor (T : TwoEI_structure) -> sig include module type of FourIdxStorage diff --git a/Basis/TwoElectronIntegralsSeparable.ml b/Basis/TwoElectronIntegralsSeparable.ml deleted file mode 100644 index a3b9156..0000000 --- a/Basis/TwoElectronIntegralsSeparable.ml +++ /dev/null @@ -1,330 +0,0 @@ -(** Two electron integral functor for operators that are separable among %{ $(x,y,z)$ %}. - It is parameterized by the [zero_m] function. -*) - -open Constants -let cutoff = integrals_cutoff - -module Bs = Basis -module Cs = ContractedShell -module Csp = ContractedShellPair -module Cspc = ContractedShellPairCouple -module Fis = FourIdxStorage - - -include FourIdxStorage - -(** Exponent of the geminal *) -let expo_s = 1.0 - -(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) -let coef_g = - [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |] - -let expo_sg_inv = - Array.map (fun x -> 1. /. (x *. expo_s *. expo_s)) - [| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |] - - - -(* - -Fit of 1/r: - -let coef_g = [| - 841.88478132 ; 70.590185207 ; 18.3616020768 ; 7.2608642093 ; - 3.57483416444 ; 2.01376031082 ; 1.24216542801 ; 0.81754348620 ; - 0.564546514023 ; 0.404228610699 ; 0.297458536575 ; 0.223321219537 ; - 0.169933732064 ; 0.130190978230 ; 0.099652303426 ; 0.075428246546 ; - 0.0555635614051 ; 0.0386791283055 ; 0.0237550435652 ; 0.010006278387 ; - |] - -let expo_sg_inv = - Array.map (fun x -> 1. /. (x *. expo_s *. expo_s)) - [| 84135.654509 ; 2971.58727634 ; 474.716025959 ; 130.676724560 ; - 47.3938388887 ; 20.2078651631 ; 9.5411021938 ; 4.8109546955 ; - 2.52795733067 ; 1.35894103210 ; 0.73586710268 ; 0.39557629706 ; - 0.20785895177 ; 0.104809693858 ; 0.049485682527 ; 0.021099788990 ; - 0.007652472186 ; 0.0021065225215 ; 0.0003365204879 ; 0.0000118855674 |] -*) - - -let class_of_contracted_shell_pair_couple shell_pair_couple = - F12RR.contracted_class_shell_pair_couple - expo_sg_inv coef_g shell_pair_couple - - - - -module Zero_m = struct - let name = "F12" -end - - - -let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs = - List.map (fun pair -> - match Cspc.make ~cutoff pair pair with - | Some cspc -> - let cls = class_of_contracted_shell_pair_couple cspc in - (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) - (* TODO \sum_k |coef_k * integral_k| *) - | None -> (pair, -1.) - ) shell_pairs - |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) - |> List.map fst - - -(* TODO - let filter_contracted_shell_pair_couples - ?(cutoff=integrals_cutoff) shell_pair_couples = - List.map (fun pair -> - let cls = - class_of_contracted_shell_pairs pair pair - in - (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) - ) shell_pairs - |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) - |> List.map fst - -*) - - -let store_class ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls = - let to_powers x = - let open Zkey in - match to_powers x with - | Three x -> x - | _ -> assert false - in - - let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple - and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple - in - - Array.iteri (fun i_c powers_i -> - let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in - let xi = to_powers powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in - let xj = to_powers powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in - let xk = to_powers powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in - let xl = to_powers powers_l in - let key = Zkey.of_powers_twelve xi xj xk xl in - let value = Zmap.find cls key in - set_chem data i_c j_c k_c l_c value - ) (Cs.zkey_array (Csp.shell_b shell_q)) - ) (Cs.zkey_array (Csp.shell_a shell_q)) - ) (Cs.zkey_array (Csp.shell_b shell_p)) - ) (Cs.zkey_array (Csp.shell_a shell_p)) - - - - - - -let of_basis_serial basis = - - let n = Bs.size basis - and shell = Bs.contracted_shells basis - in - - let eri_array = - Fis.create ~size:n `Dense -(* - Fis.create ~size:n `Sparse -*) - in - - let t0 = Unix.gettimeofday () in - - let shell_pairs = - Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs ~cutoff - in - - Printf.printf "%d significant shell pairs computed in %f seconds\n" - (List.length shell_pairs) (Unix.gettimeofday () -. t0); - - - let t0 = Unix.gettimeofday () in - let ishell = ref 0 in - - List.iter (fun shell_p -> - let () = - if (Cs.index (Csp.shell_a shell_p) > !ishell) then - (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) - in - - let sp = - Csp.shell_pairs shell_p - in - - try - List.iter (fun shell_q -> - let () = - if Cs.index (Csp.shell_a shell_q) > - Cs.index (Csp.shell_a shell_p) then - raise Exit - in - let sq = Csp.shell_pairs shell_q in - let cspc = - if Array.length sp < Array.length sq then - Cspc.make ~cutoff shell_p shell_q - else - Cspc.make ~cutoff shell_q shell_p - in - - match cspc with - | Some cspc -> - let cls = - class_of_contracted_shell_pair_couple cspc - in - store_class ~cutoff eri_array cspc cls - | None -> () - ) shell_pairs - with Exit -> () - ) shell_pairs ; - Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0); - eri_array - - - - - - - -(* Parallel functions *) - - - -let of_basis_parallel basis = - - let n = Bs.size basis - and shell = Bs.contracted_shells basis - in - - let store_class_parallel - ?(cutoff=integrals_cutoff) contracted_shell_pair_couple cls = - let to_powers x = - let open Zkey in - match to_powers x with - | Three x -> x - | _ -> assert false - in - - let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple - and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple - in - - let result = ref [] in - Array.iteri (fun i_c powers_i -> - let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in - let xi = to_powers powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in - let xj = to_powers powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in - let xk = to_powers powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in - let xl = to_powers powers_l in - let key = Zkey.of_powers_twelve xi xj xk xl in - let value = Zmap.find cls key in - result := (i_c, j_c, k_c, l_c, value) :: !result - ) (Cs.zkey_array (Csp.shell_b shell_q)) - ) (Cs.zkey_array (Csp.shell_a shell_q)) - ) (Cs.zkey_array (Csp.shell_b shell_p)) - ) (Cs.zkey_array (Csp.shell_a shell_p)); - !result - in - - - - let t0 = Unix.gettimeofday () in - - let shell_pairs = - Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs ~cutoff - in - - if Parallel.master then - Printf.printf "%d significant shell pairs computed in %f seconds\n" - (List.length shell_pairs) (Unix.gettimeofday () -. t0); - - - let t0 = Unix.gettimeofday () in - let ishell = ref max_int in - - let input_stream = Stream.of_list (List.rev shell_pairs) in - - let f shell_p = - let () = - if Parallel.rank < 2 && Cs.index (Csp.shell_a shell_p) < !ishell then - (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) - in - - let sp = - Csp.shell_pairs shell_p - in - - let result = ref [] in - try - List.iter (fun shell_q -> - let () = - if Cs.index (Csp.shell_a shell_q) > - Cs.index (Csp.shell_a shell_p) then - raise Exit - in - let sq = Csp.shell_pairs shell_q in - let cspc = - if Array.length sp < Array.length sq then - Cspc.make ~cutoff shell_p shell_q - else - Cspc.make ~cutoff shell_q shell_p - in - - match cspc with - | Some cspc -> - let cls = - class_of_contracted_shell_pair_couple cspc - in - result := (store_class_parallel ~cutoff cspc cls) :: !result; - | None -> () - ) shell_pairs; - raise Exit - with Exit -> List.concat !result |> Array.of_list - in - - let eri_array = - if Parallel.master then - Fis.create ~size:n `Dense - else - Fis.create ~size:0 `Dense - in - Farm.run ~ordered:true ~f input_stream - |> Stream.iter (fun l -> - Array.iter (fun (i_c,j_c,k_c,l_c,value) -> - set_chem eri_array i_c j_c k_c l_c value) l); - - if Parallel.master then - Printf.printf - "Computed %s Integrals in parallel in %f seconds\n%!" Zero_m.name (Unix.gettimeofday () -. t0); - Parallel.broadcast (lazy eri_array) - - - -let of_basis = - match Parallel.size with - | 1 -> of_basis_serial - | _ -> of_basis_parallel - - - - - diff --git a/Simulation.ml b/Simulation.ml index 81182b7..1b74c2f 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -4,7 +4,7 @@ type t = { nuclei : Nuclei.t; basis : Basis.t; ao_basis : AOBasis.t; - f12 : F12.f12_factor option; + f12 : F12factor.t option; nuclear_repulsion : float; } diff --git a/Simulation.mli b/Simulation.mli index 08c51ed..a7cb321 100644 --- a/Simulation.mli +++ b/Simulation.mli @@ -18,18 +18,18 @@ val ao_basis : t -> AOBasis.t val nuclear_repulsion : t -> float (** Nuclear repulsion energy *) -val f12 : t -> F12.f12_factor option +val f12 : t -> F12factor.t option (** f12 correlation factor *) (** {1 Creation} *) val make : ?cartesian:bool -> - ?multiplicity:int -> ?charge:int -> ?f12:F12.f12_factor -> + ?multiplicity:int -> ?charge:int -> ?f12:F12factor.t -> nuclei:Nuclei.t -> Basis.t -> t val of_filenames : ?cartesian:bool -> - ?multiplicity:int -> ?charge:int -> ?f12:F12.f12_factor -> + ?multiplicity:int -> ?charge:int -> ?f12:F12factor.t -> nuclei:string -> ?aux_basis_filenames:string list -> string -> t diff --git a/Utils/Util.ml b/Utils/Util.ml index 8601aea..58e00e9 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -32,6 +32,8 @@ external leadz : int64 -> int32 = "leadz_bytecode" "leadz" [@@unboxed] [@@noalloc] (** bsf instruction *) +external vfork : unit -> int = "unix_vfork" "unix_vfork" + let leadz i = leadz i |> Int32.to_int diff --git a/Utils/Util.mli b/Utils/Util.mli index 53212d7..7f083ad 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -16,6 +16,8 @@ external gamma_float : float -> float = "gamma_float_bytecode" "gamma_float" [@@unboxed] [@@noalloc] (** Gamma function [gamma] from [libm] *) +external vfork : unit -> int = "unix_vfork" "unix_vfork" + (* external popcnt : int64 -> int32 = "popcnt_bytecode" "popcnt" [@@unboxed] [@@noalloc] diff --git a/Utils/math_functions.c b/Utils/math_functions.c index 6b433ef..ba4a05a 100644 --- a/Utils/math_functions.c +++ b/Utils/math_functions.c @@ -72,3 +72,14 @@ CAMLprim value leadz_bytecode(value i) return copy_int32(__builtin_clzll (i)); } + + +#include + +CAMLprim value unix_vfork(value unit) +{ + int ret; + ret = vfork(); + return Val_int(ret); +} + diff --git a/run_fci_f12.ml b/run_fci_f12.ml index e7104f2..91e03c9 100644 --- a/run_fci_f12.ml +++ b/run_fci_f12.ml @@ -26,9 +26,11 @@ let () = arg=With_arg ""; doc="Total charge of the molecule. Default is 0"; } ; +(* { short='e' ; long="expo" ; opt=Optional; arg=With_arg ""; doc="Exponent of the Gaussian geminal"; } ; +*) { short='s' ; long="state" ; opt=Optional; arg=With_arg ""; @@ -49,11 +51,13 @@ let () = | None -> 0 in +(* let expo = match Command_line.get "expo" with | Some x -> float_of_string x | None -> 1.0 in +*) let state = match Command_line.get "state" with @@ -72,8 +76,13 @@ let () = else Printing.ppf_dev_null in +(* let f12 = - F12.gaussian_geminal expo + F12factor.gaussian_geminal expo + in +*) + let f12 = + F12factor.gaussian_geminal 1.0 in let simulation = diff --git a/run_integrals.ml b/run_integrals.ml index 98dd70e..05df36f 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -27,7 +27,7 @@ let run ~out = | Some x -> x in - let f12 = F12.gaussian_geminal 1.0 in + let f12 = F12factor.gaussian_geminal 1.0 in let s = Simulation.of_filenames ~nuclei:nuclei_file ~f12 basis_file in