diff --git a/Basis/AOBasis.ml b/Basis/AOBasis.ml index 4d0fb89..3779f91 100644 --- a/Basis/AOBasis.ml +++ b/Basis/AOBasis.ml @@ -24,7 +24,7 @@ let cartesian t = t.cartesian -let make ~cartesian ~basis nuclei = +let make ~cartesian ~basis ?f12 nuclei = let overlap = lazy ( Overlap.of_basis basis @@ -46,9 +46,15 @@ let make ~cartesian ~basis nuclei = ERI.of_basis basis ) in - let f12_ints = lazy ( - F12.of_basis basis - ) 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 { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; f12_ints ; cartesian ; diff --git a/Basis/AOBasis.mli b/Basis/AOBasis.mli index 8694136..f92d421 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 -> Nuclei.t -> t +val make : cartesian:bool -> basis:Basis.t -> ?f12:F12.f12_factor -> Nuclei.t -> t (** Creates the data structure for atomic orbitals from a {Basis.t} and the molecular geometry {Nuclei.t} *) diff --git a/Basis/F12.ml b/Basis/F12.ml index cc2cec8..22ffaa2 100644 --- a/Basis/F12.ml +++ b/Basis/F12.ml @@ -14,89 +14,86 @@ module Fis = FourIdxStorage include FourIdxStorage -(* -type gaussian_geminal = +type f12_factor = { - expo_s : float ; - coef_g : - *) -(** Exponent of the geminal *) -let expo_s = 1.0 + 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 + } -(** Slater geminal *) - -(** 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 |] +(* exp (-expo_s r) *) +let 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 * Slater *) -(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) -let coef_g = +(** 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 |] - -let expo_sg_inv = - Array.map (fun x -> 1. /. (x *. expo_s *. expo_s)) + and expo_sg = [| 0.1824 ; 0.7118; 2.252 ; 6.474 ; 19.66 ; 77.92 |] -*) - -(* -(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) -let coef_g = -[| - -3.4793465193721626604883567779324948787689208984375 ; - -0.00571703486454788484955047422886309504974633455276489257812 ; - 4.14878218728681513738365538301877677440643310546875 ; - 0.202874298181392742623785352407139725983142852783203125 ; - 0.0819187742387294803858566183407674543559551239013671875 ; - 0.04225945671351955673644695821167260874062776565551757812 ; -|] - -let expo_sg_inv = - Array.map (fun x -> 1. /. (x *. expo_s *. expo_s)) -[| - 0.63172472556807146570889699432882480323314666748046875; - 26.3759196683467962429858744144439697265625; - 0.63172102793029016876147352377302013337612152099609375; - 7.08429025944207335641067402320913970470428466796875; - 42.4442841447001910637482069432735443115234375; - 391.44036073596890901171718724071979522705078125 ; -|] -*) + in make_gaussian_corr_factor expo_s coef_g expo_sg -(* - -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 |] -*) +(* exp (-expo_s r) *) +let 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 -let class_of_contracted_shell_pair_couple shell_pair_couple = - F12RR.contracted_class_shell_pair_couple - expo_sg_inv coef_g shell_pair_couple @@ -105,13 +102,16 @@ module Zero_m = struct let name = "F12" end +let class_of_contracted_shell_pair_couple f12 shell_pair_couple = + F12RR.contracted_class_shell_pair_couple + f12.expo_sg_inv f12.coef_g shell_pair_couple -let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs = +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 cspc in + 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.) @@ -152,7 +152,7 @@ let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_coup Overlap.of_basis basis |> Overlap.matrix in - let lambda_inv = 1. /. expo_s 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 @@ -183,7 +183,7 @@ let store_class basis ?(cutoff=integrals_cutoff) data contracted_shell_pair_coup -let of_basis_serial basis = +let of_basis_serial f12 basis = let n = Bs.size basis and shell = Bs.contracted_shells basis @@ -200,7 +200,7 @@ let of_basis_serial basis = let shell_pairs = Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs ~cutoff + |> filter_contracted_shell_pairs f12 ~cutoff in Printf.printf "%d significant shell pairs computed in %f seconds\n" @@ -238,7 +238,7 @@ let of_basis_serial basis = match cspc with | Some cspc -> let cls = - class_of_contracted_shell_pair_couple cspc + class_of_contracted_shell_pair_couple f12 cspc in store_class basis ~cutoff eri_array cspc cls | None -> () @@ -258,7 +258,7 @@ let of_basis_serial basis = -let of_basis_parallel basis = +let of_basis_parallel f12 basis = let n = Bs.size basis and shell = Bs.contracted_shells basis @@ -306,7 +306,7 @@ let of_basis_parallel basis = let shell_pairs = Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs ~cutoff + |> filter_contracted_shell_pairs f12 ~cutoff in if Parallel.master then @@ -348,7 +348,7 @@ let of_basis_parallel basis = match cspc with | Some cspc -> let cls = - class_of_contracted_shell_pair_couple cspc + class_of_contracted_shell_pair_couple f12 cspc in result := (store_class_parallel ~cutoff cspc cls) :: !result; | None -> () @@ -369,7 +369,7 @@ let of_basis_parallel basis = Overlap.of_basis basis |> Overlap.matrix in - let lambda_inv = 1. /. expo_s in + let lambda_inv = 1. /. f12.expo_s in *) Farm.run ~ordered:false ~f input_stream |> Stream.iter (fun l -> @@ -377,7 +377,7 @@ let of_basis_parallel basis = (* lambda_inv *. (s.{i_c,j_c} *. s.{k_c,l_c} -. value) *) - value + value |> set_chem eri_array i_c j_c k_c l_c) l); if Parallel.master then diff --git a/CI/F12CI.ml b/CI/F12CI.ml index 781ce40..772cb49 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -2,7 +2,6 @@ open Lacaml.D type t = { - gamma : float; mo_basis : MOBasis.t ; aux_basis : MOBasis.t ; det_space : DeterminantSpace.t ; @@ -55,12 +54,6 @@ let f_ij mo_basis ki kj = in CIMatrixElement.make integrals ki kj |> List.hd - (* - match Determinant.degrees ki kj with - | (2,0) - | (0,2) -> 0.5 *. gamma *. integral - | _ -> gamma *. integral - *) let is_internal det_space = @@ -97,7 +90,7 @@ let is_internal det_space = Bitstring.logand aux_mask beta |> Bitstring.is_zero -let dressing_vector gamma ~frozen_core aux_basis f12_amplitudes ci = +let dressing_vector ~frozen_core aux_basis f12_amplitudes ci = (* let i_o1_alfa = h_ij aux_basis in @@ -215,8 +208,7 @@ let dressing_vector gamma ~frozen_core aux_basis f12_amplitudes ci = let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () = - let gamma = -0.5 /. F12.expo_s in - + let f12 = Util.of_some @@ Simulation.f12 simulation in let mo_num = MOBasis.size mo_basis in (* Add auxiliary basis set *) @@ -232,7 +224,7 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen general_basis ; GeneralBasis.read aux_basis_filename ] |> Basis.of_nuclei_and_general_basis nuclei - |> Simulation.make ~charge ~multiplicity ~nuclei + |> Simulation.make ~f12 ~charge ~multiplicity ~nuclei in let aux_basis = @@ -257,11 +249,13 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen *) ignore @@ MOBasis.f12_ints aux_basis; + let two_gamma_inv = -2. *. f12.F12.expo_s in + let f = fun ki kj -> if ki <> kj then (f_ij aux_basis ki kj) else - 1./. gamma +. (f_ij aux_basis ki kj) + two_gamma_inv +. (f_ij aux_basis ki kj) in let m_F = CI.create_matrix_spin f det_space @@ -287,7 +281,7 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen let rec iteration ?(state=1) psi = let delta = - dressing_vector gamma ~frozen_core aux_basis (f12_amplitudes psi) ci + dressing_vector ~frozen_core aux_basis (f12_amplitudes psi) ci in let f = 1.0 /. psi.{1,1} in @@ -357,6 +351,6 @@ TODO SINGLE STATE HERE ) in - { mo_basis ; aux_basis ; det_space ; ci ; eigensystem ; gamma } + { mo_basis ; aux_basis ; det_space ; ci ; eigensystem } diff --git a/Simulation.ml b/Simulation.ml index 6a3f91f..81182b7 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -4,6 +4,7 @@ type t = { nuclei : Nuclei.t; basis : Basis.t; ao_basis : AOBasis.t; + f12 : F12.f12_factor option; nuclear_repulsion : float; } @@ -13,10 +14,12 @@ let electrons t = t.electrons let basis t = t.basis let ao_basis t = t.ao_basis let nuclear_repulsion t = t.nuclear_repulsion +let f12 t = t.f12 let make ?cartesian:(cartesian=false) ?multiplicity:(multiplicity=1) ?charge:(charge=0) + ?f12 ~nuclei basis = @@ -34,7 +37,7 @@ let make ?cartesian:(cartesian=false) in let ao_basis = - AOBasis.make ~basis ~cartesian nuclei + AOBasis.make ?f12 ~basis ~cartesian nuclei in let nuclear_repulsion = @@ -42,18 +45,18 @@ let make ?cartesian:(cartesian=false) in { - charge ; basis ; nuclei ; electrons ; ao_basis ; + charge ; basis ; nuclei ; electrons ; ao_basis ; f12 ; nuclear_repulsion ; } -let of_filenames ?(cartesian=false) ?(multiplicity=1) ?(charge=0) ~nuclei ?(aux_basis_filenames=[]) basis_filename = +let of_filenames ?(cartesian=false) ?(multiplicity=1) ?(charge=0) ?f12 ~nuclei ?(aux_basis_filenames=[]) basis_filename = let nuclei = Nuclei.of_filename nuclei in let basis = Basis.of_nuclei_and_basis_filenames ~nuclei (basis_filename :: aux_basis_filenames) in - lazy (make ~cartesian ~charge ~multiplicity ~nuclei basis) + lazy (make ?f12 ~cartesian ~charge ~multiplicity ~nuclei basis) |> Parallel.broadcast diff --git a/Simulation.mli b/Simulation.mli index 947c4ba..08c51ed 100644 --- a/Simulation.mli +++ b/Simulation.mli @@ -18,14 +18,18 @@ val ao_basis : t -> AOBasis.t val nuclear_repulsion : t -> float (** Nuclear repulsion energy *) +val f12 : t -> F12.f12_factor option +(** f12 correlation factor *) + (** {1 Creation} *) val make : ?cartesian:bool -> - ?multiplicity:int -> ?charge:int -> nuclei:Nuclei.t -> Basis.t -> t + ?multiplicity:int -> ?charge:int -> ?f12:F12.f12_factor -> + nuclei:Nuclei.t -> Basis.t -> t val of_filenames : ?cartesian:bool -> - ?multiplicity:int -> ?charge:int -> nuclei:string -> - ?aux_basis_filenames:string list -> string -> t + ?multiplicity:int -> ?charge:int -> ?f12:F12.f12_factor -> + nuclei:string -> ?aux_basis_filenames:string list -> string -> t diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index 490e3ae..b08f1bc 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -358,7 +358,7 @@ let four_index_transform coef source = 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); + if Parallel.master then (incr n ; Printf.eprintf "\r%d / %d%!" !n mo_num); Array.iter (fun (alpha, beta, gamma, delta, x) -> set_chem destination alpha beta gamma delta x) l); diff --git a/run_fci_f12.ml b/run_fci_f12.ml index aed25c9..c65a598 100644 --- a/run_fci_f12.ml +++ b/run_fci_f12.ml @@ -23,6 +23,10 @@ let () = { short='c' ; long="charge" ; opt=Optional; 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"; } ; ] end; @@ -39,6 +43,12 @@ let () = | None -> 0 in + let expo = + match Command_line.get "expo" with + | Some x -> float_of_string x + | None -> 1.0 + in + let multiplicity = match Command_line.get "multiplicity" with | Some x -> int_of_string x @@ -50,8 +60,12 @@ let () = else Printing.ppf_dev_null in + let f12 = + F12.gaussian_geminal expo + in + let simulation = - Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file + Simulation.of_filenames ~f12 ~charge ~multiplicity ~nuclei:nuclei_file basis_file in let hf = HartreeFock.make simulation in @@ -64,7 +78,9 @@ let () = let fcif12 = F12CI.make ~simulation ~frozen_core:true ~mo_basis ~aux_basis_filename () in + let ci = F12CI.ci fcif12 in Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion simulation); + let _, e_cif12 = F12CI.eigensystem fcif12 in Format.fprintf ppf "FCI-F12 energy : %20.16f@." (e_cif12.{1} +. Simulation.nuclear_repulsion simulation);