From 08df6eddec4c1916e3bdb365849b3343278ab8c6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 2 Feb 2018 21:49:07 +0100 Subject: [PATCH] Recursive Boys function + removed zerom cache --- Basis/TwoElectronRR.ml | 16 +--- Basis/TwoElectronRRVectorized.ml | 29 +------ Utils/Util.ml | 141 ++++++++++++++++--------------- 3 files changed, 74 insertions(+), 112 deletions(-) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index defbd5b..8e0a7dc 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -325,9 +325,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (* Compute all integrals in the shell for each pair of significant shell pairs *) - let zero_m_cache = - Hashtbl.create 129 - in for ab=0 to (Array.length shell_p - 1) do let cab = shell_p.(ab).Shell_pair.coef in let b = shell_p.(ab).Shell_pair.j in @@ -352,18 +349,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let zero_m_array = - let key = (maxm, expo_pq_inv, norm_pq_sq) in - try - let result = - Hashtbl.find zero_m_cache key - in - result - with - | Not_found -> - let result = - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq - in - (Hashtbl.add zero_m_cache key result ; result) + zero_m ~maxm ~expo_pq_inv ~norm_pq_sq in begin match Contracted_shell.(totAngMom shell_a, totAngMom shell_b, diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 48e62af..6123a8c 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -306,9 +306,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.make (Array.length class_indices) 0.; in - let zero_m_cache = - Hashtbl.create 129 - in (* Compute all integrals in the shell for each pair of significant shell pairs *) begin @@ -338,18 +335,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let zero_m_array = - let key = (0, expo_pq_inv, norm_pq_sq) in - try - let result = - Hashtbl.find zero_m_cache key - in - result - with - | Not_found -> - let result = - zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq - in - (Hashtbl.add zero_m_cache key result ; result) + zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq in accu +. coef_prod *. zero_m_array.(0) @@ -377,18 +363,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let zero_m_array = - let key = (maxm, expo_pq_inv, norm_pq_sq) in - try - let result = - Hashtbl.find zero_m_cache key - in - result - with - | Not_found -> - let result = - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq - in - (Hashtbl.add zero_m_cache key result ; result) + zero_m ~maxm ~expo_pq_inv ~norm_pq_sq in let d = shell_cd.Shell_pair.j in diff --git a/Utils/Util.ml b/Utils/Util.ml index 258b028..330820b 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -10,86 +10,53 @@ let factmax = 150 -(*reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten + +(* 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 + + reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten (New Algorithm handbook in C language) (Gijyutsu hyouron sha, Tokyo, 1991) p.227 [in Japanese] *) -(* Incomplete gamma function - 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 *) - -let rec p_gamma a x loggamma_a = - if x >= 1. +. a then 1. -. q_gamma a x loggamma_a - else if x = 0. then 0. - else - let rec pg_loop prev res term k = - if k > 1000. then failwith "p_gamma did not converge." - else if prev = res then res - else - let term = term *. x /. (a +. k) in - pg_loop res (res +. term) term (k +. 1.) - in - let r0 = exp (a *. log x -. x -. loggamma_a) /. a in - pg_loop min_float r0 r0 1. - - and q_gamma a x loggamma_a = - if x < 1. +. a then 1. -. p_gamma a x loggamma_a +let incomplete_gamma ~alpha x = + let rec p_gamma a x loggamma_a = + if x >= 1. +. a then 1. -. q_gamma a x loggamma_a + else if x = 0. then 0. else - let rec qg_loop prev res la lb w k = - if k > 1000. then failwith "q_gamma did not converge." + let rec pg_loop prev res term k = + if k > 1000. then failwith "p_gamma did not converge." else if prev = res then res else - let k_inv = 1. /. k in - let la, lb = - lb, ((k -. 1. -. a) *. (lb -. la) +. (k +. x) *. lb) *. k_inv - in - let w = w *. (k -. 1. -. a) *. k_inv in - let prev, res = res, res +. w /. (la *. lb) in - qg_loop prev res la lb w (k +. 1.) + let term = term *. x /. (a +. k) in + pg_loop res (res +. term) term (k +. 1.) in - let w = exp (a *. log x -. x -. loggamma_a) in - let lb = (1. +. x -. a) in - qg_loop min_float (w /. lb) 1. lb w 2.0 + let r0 = exp (a *. log x -. x -. loggamma_a) /. a in + pg_loop min_float r0 r0 1. - -(** Generalized Boys function. Uses GSL's incomplete Gamma function. - maxm : Maximum total angular momentum -*) -let rec boys_function ~maxm t = - match maxm with - | 0 -> - begin - if t = 0. then [| 1. |] else - let sq_t = sqrt t in - [| (sq_pi_over_two /. sq_t) *. erf_float sq_t |] - end - | _ -> - begin - if t = 0. then - Array.init (maxm+1) (fun m -> 1. /. float_of_int (m+m+1)) + and q_gamma a x loggamma_a = + if x < 1. +. a then 1. -. p_gamma a x loggamma_a else - let incomplete_gamma ~alpha x = - let gf = gamma_float alpha in - gf *. p_gamma alpha x (log gf) - in - let t_inv = - 1. /. t - in - let factor = - Array.make (maxm+1) (0.5, sqrt t_inv); - in - for i=1 to maxm - do - let (dm, f) = factor.(i-1) in - factor.(i) <- (dm +. 1., f *. t_inv); - done; - Array.map (fun (dm, f) -> - if dm = 0.5 then - (boys_function ~maxm:0 t).(0) + let rec qg_loop prev res la lb w k = + if k > 1000. then failwith "q_gamma did not converge." + else if prev = res then res else - (incomplete_gamma dm t ) *. 0.5 *. f - ) factor - end + let k_inv = 1. /. k in + let la, lb = + lb, ((k -. 1. -. a) *. (lb -. la) +. (k +. x) *. lb) *. k_inv + in + let w = w *. (k -. 1. -. a) *. k_inv in + let prev, res = res, res +. w /. (la *. lb) in + qg_loop prev res la lb w (k +. 1.) + in + let w = exp (a *. log x -. x -. loggamma_a) in + let lb = (1. +. x -. a) in + qg_loop min_float (w /. lb) 1. lb w 2.0 + in + let gf = gamma_float alpha in + gf *. p_gamma alpha x (log gf) + + @@ -139,3 +106,37 @@ let chop f g = +(** Generalized Boys function. + maxm : Maximum total angular momentum +*) +let boys_function ~maxm t = + match maxm with + | 0 -> + begin + if t = 0. then [| 1. |] else + let sq_t = sqrt t in + [| (sq_pi_over_two /. sq_t) *. erf_float sq_t |] + end + | _ -> + begin + let result = + Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1)) + in + if t <> 0. then + begin + let fmax = + let t_inv = sqrt (1. /. t) in + let n = float_of_int maxm in + let dm = 0.5 +. n in + let f = (pow t_inv (maxm+maxm+1) ) in + (incomplete_gamma dm t) *. 0.5 *. f + in + let emt = exp (-. t) in + result.(maxm) <- fmax; + for n=maxm-1 downto 0 do + result.(n) <- ( (t+.t) *. result.(n+1) +. emt) *. result.(n) + done + end; + result + end +