From 4e8efef26441572aecca40d2bd9542182941fa41 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 17 Feb 2020 19:45:53 +0100 Subject: [PATCH 1/2] Added ScreenedERI --- Basis/AOBasis.ml | 22 +++++--- Basis/F12.ml | 2 +- Basis/F12RR.ml | 11 ++-- Basis/GaussianOperator.ml | 18 +++--- Basis/KinInt.ml | 4 +- Basis/OneElectronRR.ml | 2 +- Basis/Overlap_primitives.ml | 2 +- Basis/ScreenedERI.ml | 95 ++++++++++++++++++++++++++++++++ Basis/TwoElectronIntegrals.ml | 8 +-- Basis/TwoElectronRR.ml | 25 +++++---- Basis/TwoElectronRRVectorized.ml | 29 +++++----- Basis/Zero_m_parameters.ml | 6 +- Utils/Util.ml | 48 +++++++++++++--- Utils/Util.mli | 10 +++- 14 files changed, 216 insertions(+), 66 deletions(-) create mode 100644 Basis/ScreenedERI.ml diff --git a/Basis/AOBasis.ml b/Basis/AOBasis.ml index c4700cb..83b678b 100644 --- a/Basis/AOBasis.ml +++ b/Basis/AOBasis.ml @@ -3,13 +3,14 @@ open Util type t = { - basis : Basis.t ; - overlap : Overlap.t lazy_t; - ortho : Orthonormalization.t lazy_t; - eN_ints : NucInt.t lazy_t; - kin_ints : KinInt.t lazy_t; - ee_ints : ERI.t lazy_t; - f12_ints : F12.t lazy_t; + basis : Basis.t ; + overlap : Overlap.t lazy_t; + ortho : Orthonormalization.t lazy_t; + eN_ints : NucInt.t lazy_t; + kin_ints : KinInt.t lazy_t; + ee_ints : ERI.t lazy_t; + f12_ints : F12.t lazy_t; + f12_over_r12_ints : ScreenedERI.t lazy_t; cartesian : bool; } @@ -20,6 +21,7 @@ let eN_ints t = Lazy.force t.eN_ints let kin_ints t = Lazy.force t.kin_ints let ee_ints t = Lazy.force t.ee_ints let f12_ints t = Lazy.force t.f12_ints +let f12_over_r12_ints t = Lazy.force t.f12_over_r12_ints let cartesian t = t.cartesian @@ -50,7 +52,11 @@ let make ~cartesian ~basis ?f12 nuclei = F12.of_basis basis ) in - { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; f12_ints ; + let f12_over_r12_ints = lazy ( + ScreenedERI.of_basis basis + ) in + + { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; f12_ints ; f12_over_r12_ints ; cartesian ; } diff --git a/Basis/F12.ml b/Basis/F12.ml index c24b64b..1cf2f7b 100644 --- a/Basis/F12.ml +++ b/Basis/F12.ml @@ -14,7 +14,7 @@ module T = struct let class_of_contracted_shell_pair_couple shell_pair_couple = let g = f12_factor.F12factor.gaussian in F12RR.contracted_class_shell_pair_couple - g.GaussianOperator.expo_sg_inv + g.GaussianOperator.expo_g_inv g.GaussianOperator.coef_g shell_pair_couple diff --git a/Basis/F12RR.ml b/Basis/F12RR.ml index 90fc56a..ce0f353 100644 --- a/Basis/F12RR.ml +++ b/Basis/F12RR.ml @@ -12,7 +12,6 @@ module Po = Powers module Psp = PrimitiveShellPair module Pspc = PrimitiveShellPairCouple module Ps = PrimitiveShell -module Zp = Zero_m_parameters let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff @@ -83,7 +82,7 @@ let rec hvrr angMom_a angMom_b angMom_c angMom_d let v1 = vrr (angMom_a-1) 0 in - let a = float_of_int (angMom_a-1) in + let a = float_of_int_fast (angMom_a-1) in let v2 = vrr (angMom_a-2) 0 in @@ -92,7 +91,7 @@ let rec hvrr angMom_a angMom_b angMom_c angMom_d let v1 = vrr angMom_a 0 in - let a = float_of_int angMom_a in + let a = float_of_int_fast angMom_a in let v2 = vrr (angMom_a-1) 0 in @@ -102,7 +101,7 @@ let rec hvrr angMom_a angMom_b angMom_c angMom_d let v1 = vrr 0 (angMom_c-1) in - let b = float_of_int (angMom_c-1) in + let b = float_of_int_fast (angMom_c-1) in let v3 = vrr 0 (angMom_c-2) in @@ -112,8 +111,8 @@ let rec hvrr angMom_a angMom_b angMom_c angMom_d let v1 = vrr angMom_a (angMom_c-1) in - let a = float_of_int angMom_a in - let b = float_of_int (angMom_c-1) in + let a = float_of_int_fast angMom_a in + let b = float_of_int_fast (angMom_c-1) in let v2 = vrr (angMom_a-1) (angMom_c-1) in diff --git a/Basis/GaussianOperator.ml b/Basis/GaussianOperator.ml index 01a1242..717525f 100644 --- a/Basis/GaussianOperator.ml +++ b/Basis/GaussianOperator.ml @@ -4,16 +4,16 @@ open Constants type t = { - coef_g : float array; - expo_sg : float array; - expo_sg_inv : float array; + coef_g : float array; + expo_g : float array; + expo_g_inv : float array; } -let make coef_g expo_sg = - let expo_sg_inv = - Array.map (fun x -> 1. /. x ) expo_sg +let make coef_g expo_g = + let expo_g_inv = + Array.map (fun x -> 1. /. x ) expo_g in - { coef_g ; expo_sg ; expo_sg_inv } + { coef_g ; expo_g ; expo_g_inv } let one_over_r = @@ -25,13 +25,13 @@ let one_over_r = 0.169933732064 ; 0.130190978230 ; 0.099652303426 ; 0.075428246546 ; 0.0555635614051 ; 0.0386791283055 ; 0.0237550435652 ; 0.010006278387 ; |] - and expo_sg = + and expo_g = [| 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 + in make coef_g expo_g diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index 8472721..1f849c1 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -90,8 +90,8 @@ let contracted_class shell_a shell_b : float Zmap.t = and s2 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB - 1) k and s3 = ov (Po.get xyz angMomA - 1) (Po.get xyz angMomB + 1) k and s4 = ov (Po.get xyz angMomA + 1) (Po.get xyz angMomB + 1) k - and a = float_of_int (Po.get xyz angMomA) - and b = float_of_int (Po.get xyz angMomB) + and a = float_of_int_fast (Po.get xyz angMomA) + and b = float_of_int_fast (Po.get xyz angMomB) in 0.5 *. a *. b *. s1 -. expo_a *. b *. s2 -. expo_b *. a *. s3 +. 2.0 *. expo_a *. expo_b *. s4 diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index 96bb9cb..8b850dd 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -58,7 +58,7 @@ let hvrr_one_e angMom_a angMom_b vrr amm in let v1 = vrr am in - let f3 = float_of_int amxyz *. expo_inv_p *. 0.5 in + let f3 = float_of_int_fast amxyz *. expo_inv_p *. 0.5 in Array.init (maxsze - angMom_a.Po.tot) (fun m -> f1 *. v1.(m) -. f2 *. v1.(m+1) +. f3 *. (v3.(m) -. expo_inv_p *. v3.(m+1)) diff --git a/Basis/Overlap_primitives.ml b/Basis/Overlap_primitives.ml index 41eb0c3..93c4065 100644 --- a/Basis/Overlap_primitives.ml +++ b/Basis/Overlap_primitives.ml @@ -10,7 +10,7 @@ let hvrr (angMom_a, angMom_b) (expo_inv_p) (center_ab, center_pa) if angMom_a > 0 then chop center_pa (fun () -> vrr (angMom_a-1)) - +. chop (0.5 *. float_of_int (angMom_a-1) *. expo_inv_p) + +. chop (0.5 *. float_of_int_fast (angMom_a-1) *. expo_inv_p) (fun () -> vrr (angMom_a-2)) else if angMom_a = 0 then 1. else 0. diff --git a/Basis/ScreenedERI.ml b/Basis/ScreenedERI.ml new file mode 100644 index 0000000..14c482e --- /dev/null +++ b/Basis/ScreenedERI.ml @@ -0,0 +1,95 @@ +(** Screened Electron-electron repulsion integrals (Yukawa potential). Required for F12/r12. *) + +open Constants +open Util + +module Csp = ContractedShellPair +module Cspc = ContractedShellPairCouple + +module T = struct + + let name = "Screened electron repulsion integrals" + + let f12_factor = F12factor.gaussian_geminal 1.0 + + open Zero_m_parameters + + let zero_m z = + let expo_pq_inv = z.expo_p_inv +. z.expo_q_inv in + assert (expo_pq_inv <> 0.); + (* + let f12_factor = + match z.f12_factor with + | Some f -> f + | None -> assert false + in + *) + + let expo_G_inv, coef_G = + f12_factor.F12factor.gaussian.GaussianOperator.expo_g_inv, + f12_factor.F12factor.gaussian.GaussianOperator.coef_g + in + + let expo_pq = 1. /. expo_pq_inv in + let maxm = z.maxm in + + let result = Array.init (maxm+1) (fun _ -> 0.) in + array_range 0 (Array.length coef_G) + |> Array.iter (fun i -> + let fG_inv = expo_pq_inv +. expo_G_inv.(i) in + let fG = 1. /. fG_inv in + let t = + if z.norm_pq_sq > integrals_cutoff then + z.norm_pq_sq *. (expo_pq -. fG) + else 0. + in + let fm = boys_function ~maxm t in + + let tmp_array = + let result = Array.init (maxm+1) (fun k -> 1.) in + array_range 1 maxm + |> Array.iter (fun k -> result.(k) <- result.(k-1) *. expo_pq *. expo_G_inv.(i)); + result + in + + let tmp_result = Array.init (maxm+1) (fun _ -> 0.) in + let rec aux accu m = function + | 0 -> tmp_result.(m) <- fm.(m) *. accu + | j -> + begin + let f = + array_range 0 m + |> Array.fold_left (fun v k -> + v +. (binom_float m k) *. tmp_array.(k) *. fm.(k)) 0. + in + tmp_result.(m) <- accu *. f; + let new_accu = -. accu *. expo_pq in + (aux [@tailcall]) new_accu (m+1) (j-1) + end + in + let f = + two_over_sq_pi *. (sqrt expo_pq) *. fG *. expo_G_inv.(i) *. exp (-.fG *. z.norm_pq_sq) + in + aux f 0 maxm; + Array.iteri (fun k v -> + result.(k) <- result.(k) +. coef_G.(i) *. v + ) tmp_result) ; + 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 = TwoElectronIntegrals.Make(T) +include M + diff --git a/Basis/TwoElectronIntegrals.ml b/Basis/TwoElectronIntegrals.ml index 7d4530d..4fe99b1 100644 --- a/Basis/TwoElectronIntegrals.ml +++ b/Basis/TwoElectronIntegrals.ml @@ -27,7 +27,7 @@ module Make(T : TwoEI_structure) = struct 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 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.) @@ -96,9 +96,9 @@ module Make(T : TwoEI_structure) = struct let eri_array = Fis.create ~size:n `Dense - (* - Fis.create ~size:n `Sparse - *) + (* + Fis.create ~size:n `Sparse + *) in let t0 = Unix.gettimeofday () in diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 2a51148..ebc73b6 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -15,6 +15,7 @@ module Zp = Zero_m_parameters let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff + exception NullQuartet @@ -105,7 +106,7 @@ let rec hvrr_two_e let amm = Po.decr xyz am in let v3 = vrr0 amm in let v1 = vrr0 am in - let f3 = (float_of_int amxyz) *. expo_p_inv *. 0.5 in + let f3 = (float_of_int_fast amxyz) *. expo_p_inv *. 0.5 in Array.iteri (fun m _ -> result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) +. f3 *. (v3.(m) +. expo_p_inv *. v3.(m+1)) ) result @@ -142,7 +143,7 @@ let rec hvrr_two_e begin let am = Po.decr xyz angMom_a in let f5 = - (float_of_int axyz) *. expo_p_inv *. expo_q_inv *. 0.5 + (float_of_int_fast axyz) *. expo_p_inv *. expo_q_inv *. 0.5 in if (abs_float f5 > cutoff) then let v5 = @@ -154,7 +155,7 @@ let rec hvrr_two_e if cmxyz > 0 then begin let f3 = - (float_of_int cmxyz) *. expo_q_inv *. 0.5 + (float_of_int_fast cmxyz) *. expo_q_inv *. 0.5 in if (abs_float f3 > cutoff) || (abs_float (f3 *. expo_q_inv) > cutoff) then @@ -206,7 +207,7 @@ let rec hvrr_two_e let result = if cmxyz < 1 then result else - let f = 0.5 *. (float_of_int cmxyz) *. expo_q_inv in + let f = 0.5 *. (float_of_int_fast cmxyz) *. expo_q_inv in if abs_float f < cutoff then 0. else let cmm = Po.decr xyz cm in let v3 = trr angMom_a cmm in @@ -226,7 +227,7 @@ let rec hvrr_two_e in let result = if axyz < 1 then result else - let f = 0.5 *. (float_of_int axyz) *. expo_q_inv in + let f = 0.5 *. (float_of_int_fast axyz) *. expo_q_inv in if abs_float f < cutoff then result else let am = Po.decr xyz angMom_a in let v2 = trr am cm in @@ -333,17 +334,16 @@ let contracted_class_shell_pair_couple ~zero_m shell_pair_couple : float Zmap.t and sp_cd = Pspc.shell_pair_q spc in - let expo_p_inv = Psp.exponent_inv sp_ab - in - let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in let center_pa = Psp.center_minus_a sp_ab in let center_qc = Psp.center_minus_a sp_cd in let norm_pq_sq = Co.dot center_pq center_pq in + let expo_p_inv = Psp.exponent_inv sp_ab in let expo_q_inv = Psp.exponent_inv sp_cd in let normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in - let zero_m_array = - zero_m Zp.{ + let zero = Zp.zero zero_m in + let zero_m_array = zero_m + { zero with maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; center_pq ; center_pa ; center_qc ; normalization ; } @@ -456,8 +456,9 @@ let contracted_class_atomic_shell_pair_couple ~zero_m atomic_shell_pair_couple : let expo_q_inv = Psp.exponent_inv sp_cd in let normalization = Psp.normalization sp_ab *. Psp.normalization sp_cd in - let zero_m_array = - zero_m Zp.{ + let zero = Zp.zero zero_m in + let zero_m_array = zero_m + { zero with maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; center_pq ; center_pa ; center_qc ; normalization ; } diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 6e5c8d4..293f45c 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -112,7 +112,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) else begin let amm = Po.decr xyz am in - let amxyz = float_of_int amxyz in + let amxyz = float_of_int_fast amxyz in let v_amm = vrr0_v amm in Array.iteri (fun l expo_inv_p_l -> let f = amxyz *. expo_p_inv.(l) *. 0.5 @@ -231,7 +231,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let p2 = if cxyz < 2 then p1 else let cmm = Po.decr xyz cm in - let fcm = (float_of_int (cxyz-1)) *. 0.5 in + let fcm = (float_of_int_fast (cxyz-1)) *. 0.5 in let f1 = Array.init nq (fun k -> let x = fcm *. expo_q_inv.(k) in @@ -356,7 +356,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | Some p2 -> p2 in for l=0 to np-1 do - let fa = (float_of_int axyz) *. expo_p_inv.(l) *. 0.5 in + let fa = (float_of_int_fast axyz) *. expo_p_inv.(l) *. 0.5 in let p2_l = p2.(l) and v_l = v.(l) in @@ -396,7 +396,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result = if cmxyz < 1 then result else begin - let f = 0.5 *. (float_of_int cmxyz) in + let f = 0.5 *. (float_of_int_fast cmxyz) in let cmm = Po.decr xyz cm in match result, trr_v angMom_a cmm with | None, None -> None @@ -473,7 +473,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result = if axyz < 1 then result else begin - let f = 0.5 *. (float_of_int axyz) in + let f = 0.5 *. (float_of_int_fast axyz) in let am = Po.decr xyz angMom_a in match result, trr_v am cm with | Some result, None -> Some result @@ -655,9 +655,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Psp.normalization sq.(i-1) in + let zero = Zp.zero zero_m in let zero_m_array = - zero_m Zp.{ - maxm=0 ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; + zero_m + {zero with + expo_p_inv ; expo_q_inv ; norm_pq_sq ; center_pq ; center_pa ; center_qc ; normalization ; } in @@ -787,12 +789,13 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let normalization = Psp.normalization shell_ab *. Psp.normalization shell_cd in - zero_m Zp.{ - maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; - center_pq = Coordinate.make Coordinate.{x ; y ; z} ; - center_pa ; center_qc = center_qc_tmp.(cd) ; - normalization ; - } + let zero = Zp.zero zero_m in + zero_m {zero with + maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ; + center_pq = Coordinate.make Coordinate.{x ; y ; z} ; + center_pa ; center_qc = center_qc_tmp.(cd) ; + normalization ; + } ) sq in (* Transpose result *) diff --git a/Basis/Zero_m_parameters.ml b/Basis/Zero_m_parameters.ml index 0d8ebd8..d19ef6c 100644 --- a/Basis/Zero_m_parameters.ml +++ b/Basis/Zero_m_parameters.ml @@ -8,10 +8,14 @@ type t = center_pq : Coordinate.t ; center_pa : Coordinate.t ; center_qc : Coordinate.t ; + f12_factor : F12factor.t option; + zero_m_func : t -> float array ; } -let zero = +let zero ?f12_factor zero_m_func = { + zero_m_func ; + f12_factor ; maxm=0 ; expo_p_inv = 0.; expo_q_inv = 0.; diff --git a/Utils/Util.ml b/Utils/Util.ml index 6744357..8482ddd 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -45,12 +45,22 @@ let () = ;; +let memo_float_of_int = + Array.init 64 float_of_int + +let float_of_int_fast i = + if Int.logand i 63 = i then + memo_float_of_int.(i) + else + float_of_int i + let factmax = 150 -(* Incomplete gamma function : Int_0^x exp(-t) t^(a-1) dt - p: 1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt - q: 1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt +(* Incomplete gamma function : + {% $\gamma(\alpha,x) = \int_0^x e^{-t} t^{\alpha-1} dt$ %} + {% $p: \frac{1}{\Gamma(\alpha)} \int_0^x e^{-t} t^{\alpha-1} dt$ %} + {% $q: \frac{1}{\Gamma(\alpha)} \int_x^\infty e^{-t} t^{\alpha-1} dt$ %} reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten (New Algorithm handbook in C language) (Gijyutsu hyouron @@ -123,7 +133,7 @@ let fact = function | i -> fact_memo.(i) -let binom n k = +let binom_func n k = (*TODO : slow function *) assert (n >= k); let rec aux n k = @@ -133,6 +143,23 @@ let binom n k = aux (n-1) (k-1) + aux (n-1) k in aux n k +let binom = + let memo = + Array.init 64 (fun n -> Array.init 64 (fun k -> + if n >= k then binom_func n k else 0) ) + in + fun n k -> + if n < 64 then + begin + assert (n >= k); + memo.(n).(k) + end + else + binom_func n k + +let binom_float n k = + binom n k + |> float_of_int_fast let rec pow a = function | 0 -> 1. @@ -157,6 +184,14 @@ let chop f g = (** Generalized Boys function. maxm : Maximum total angular momentum + {% $F_m(x) = \frac{\gamma(m+1/2,x)}{2x^{m+1/2}}$ %} + where %{ $\gamma$ %} is the incomplete gamma function. + + {% $F_0(0.) = 1$ %} + {% $F_0(t) = \frac{\sqrt{\pi}}{2\sqrt{t}} \text{erf} ( \sqrt{t} )$ %} + {% $F_m(0.) = \frac{1}{2m+1}$ %} + {% $F_m(t) = \frac{\gamma{m+1/2,t}}{2t^{m+1/2}} + {% $F_m(t) = \frac{ 2t\, F_{m+1}(t) + e^{-t} }{2m+1}$ %} *) let boys_function ~maxm t = match maxm with @@ -182,10 +217,9 @@ let boys_function ~maxm t = let dm = 0.5 +. n in let f = (pow t_inv (maxm+maxm+1) ) in match classify_float f with + | FP_normal -> (incomplete_gamma dm t) *. 0.5 *. f | FP_zero - | FP_subnormal - | FP_normal -> - (incomplete_gamma dm t) *. 0.5 *. f + | FP_subnormal -> 0. | _ -> invalid_arg "zero_m overflow" in let emt = exp (-. t) in diff --git a/Utils/Util.mli b/Utils/Util.mli index 8aedfa5..c76599a 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -49,6 +49,9 @@ val fact : int -> float val binom : int -> int -> int (** Binomial coefficient. [binom n k] {% $= C_n^k$ %}. *) +val binom_float : int -> int -> float +(** Binomial coefficient. [binom n k] {% $= C_n^k$ %}. *) + val pow : float -> int -> float (** Fast implementation of the power function for small integer powers *) @@ -57,8 +60,10 @@ val chop : float -> (unit -> float) -> float than {!Constants.epsilon}, and return [a *. f ()]. *) -val of_some : 'a option -> 'a +val float_of_int_fast : int -> float +(* Faster implementation of float_of_int for small positive ints *) +val of_some : 'a option -> 'a (** {2 Functions related to the Boys function} *) @@ -77,6 +82,9 @@ val boys_function : maxm:int -> float -> float array (** {2 Extension of the Array module} *) +val array_range : int -> int -> int array +(** [array_range first last] returns an array [| first; first+1 ; ... ; last-1 ; last |]. *) + val array_sum : float array -> float (** Returns the sum of all the elements of the array *) From 873342f8edb796884b0071ea112458eebfd34770 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 20 Feb 2020 23:27:45 +0100 Subject: [PATCH 2/2] Improved FCI program --- run_fci.ml | 51 +++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 41 insertions(+), 10 deletions(-) diff --git a/run_fci.ml b/run_fci.ml index 6b3dc4a..72c4d38 100644 --- a/run_fci.ml +++ b/run_fci.ml @@ -1,8 +1,10 @@ +open Lacaml.D + 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 Full CI calculation"; set_specs [ { short='b' ; long="basis" ; opt=Mandatory; arg=With_arg ""; @@ -19,6 +21,14 @@ let () = { short='c' ; long="charge" ; opt=Optional; arg=With_arg ""; doc="Total charge of the molecule. Default is 0"; } ; + + { short='f' ; long="frozen-core" ; opt=Optional; + arg=Without_arg ; + doc="Freeze core MOs. Default is false."; } ; + + { short='s' ; long="state" ; opt=Optional; + arg=With_arg ""; + doc="Requested state. Default is 1."; } ; ] end; @@ -29,10 +39,19 @@ let () = let charge = match Command_line.get "charge" with - | Some x -> int_of_string x + | Some x -> ( if x.[0] = 'm' then + ~- (int_of_string (String.sub x 1 (String.length x - 1))) + else + int_of_string x ) | None -> 0 in + let state = + match Command_line.get "state" with + | Some x -> int_of_string x + | None -> 1 + in + let multiplicity = match Command_line.get "multiplicity" with | Some x -> int_of_string x @@ -44,25 +63,37 @@ 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 - Format.fprintf ppf "@[%a@]@." HartreeFock.pp hf; + let hf = HartreeFock.make simulation in + if Parallel.master then + Format.fprintf ppf "@[%a@]@." HartreeFock.pp hf; - let mos = + let mo_basis = MOBasis.of_hartree_fock hf in - let space = - DeterminantSpace.fci_of_mo_basis ~frozen_core:false mos + let frozen_core = + Command_line.get_bool "frozen-core" in - let ci = CI.make space in - Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); + + let space = + DeterminantSpace.fci_of_mo_basis ~frozen_core mo_basis + in + + let ci = CI.make space ~n_states:state in + if Parallel.master then + Format.fprintf ppf "FCI energy : "; + Vec.iteri (fun i x -> if i <= state then + Format.fprintf ppf "%20.16f@; " (x +. Simulation.nuclear_repulsion simulation) ) + (CI.eigenvalues ci); + (* let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in Util.list_range 1 (DeterminantSpace.size space) |> List.iter (fun i -> Format.printf "@[%f@]@;" s2.{i,i}); *) + Format.fprintf ppf "@.";