LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2025-01-03 01:55:40 +01:00
Code Releases Activity

Fast 00

Browse Source
This commit is contained in:
Anthony Scemama 2018-02-01 19:06:29 +01:00
parent 938e9c18fb
commit c666878dcd
2 changed files with 30 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Utils/Constants.ml
View File

@ -2,6 +2,8 @@ let cutoff = 1.e-15
(** Constants *)
let pi = acos (-1.)
let sq_pi = sqrt pi
let sq_pi_over_two = sq_pi *. 0.5
let pi_inv = 1. /. pi
let two_over_sq_pi = 2. /. (sqrt pi)

48
Utils/Util.ml
View File

@ -7,26 +7,34 @@ let factmax = 150
maxm : Maximum total angular momentum
*)
let boys_function ~maxm t =
begin
if t = 0. then
Array.init (maxm+1) (fun m -> 1. /. float_of_int (m+m+1))
else
let incomplete_gamma ~alpha x =
Gsl.Sf.gamma alpha *. ( Gsl.Sf.gamma_inc_P alpha x )
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) -> (incomplete_gamma dm t ) *. 0.5 *. f) factor
end
match maxm with
| 0 ->
begin
if t = 0. then [| 1. |] else
let sq_t = sqrt t in
[| (sq_pi_over_two /. sq_t) *. Gsl.Sf.erf sq_t |]
end
| _ ->
begin
if t = 0. then
Array.init (maxm+1) (fun m -> 1. /. float_of_int (m+m+1))
else
let incomplete_gamma ~alpha x =
Gsl.Sf.gamma alpha *. ( Gsl.Sf.gamma_inc_P alpha x )
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) -> (incomplete_gamma dm t ) *. 0.5 *. f) factor
end