mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-07-25 04:07:24 +02:00
Recursive Boys function + removed zerom cache
This commit is contained in:
parent
7ee34f3c9f
commit
08df6eddec
@ -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,
|
||||
|
||||
@ -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
|
||||
|
||||
141
Utils/Util.ml
141
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
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user