LCPQ
/
QCaml
mirror of https://gitlab.com/scemama/QCaml.git
10
1
Fork 0
Code Releases Activity

Accelerated norm function

This commit is contained in:
Anthony Scemama 2018-02-01 16:09:04 +01:00
parent 1aeb890aa9
commit f3d03cd424
4 changed files with 161 additions and 115 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
Basis/Contracted_shell.ml
View File

@ -6,7 +6,8 @@ type t = {
center : Coordinate.t;
totAngMom : Angular_momentum.t;
size : int;
norm_coef : (int array -> float) array;
norm_coef : float array;
norm_coef_scale : float Zmap.t;
indice : int;
powers : Zkey.t array;
}
@ -17,6 +18,7 @@ let coef a i = a.coef.(i)
let center a = a.center
let totAngMom a = a.totAngMom
let norm_coef a i = a.norm_coef.(i)
let norm_coef_scale a = a.norm_coef_scale
let indice a = a.indice
let powers a = a.powers
@ -71,10 +73,20 @@ let compute_norm_coef expo totAngMom =
let create ~indice ~expo ~coef ~center ~totAngMom =
assert (Array.length expo = Array.length coef);
assert (Array.length expo > 0);
let norm_coef =
let norm_coef_func =
compute_norm_coef expo totAngMom
in
let powers =
Angular_momentum.zkey_array (Angular_momentum.Singlet totAngMom)
in
let norm_coef =
Array.map (fun f -> f [| Angular_momentum.to_int totAngMom ; 0 ; 0 |]) norm_coef_func
in
let norm_coef_scale =
Zmap.create 13
in
Array.iter (fun a -> Zmap.add norm_coef_scale a ((norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0))) powers;
{ indice ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
powers = Angular_momentum.zkey_array (Angular_momentum.Singlet totAngMom) }
norm_coef_scale ; powers }

41
Basis/Shell_pair.ml
View File

@ -7,7 +7,7 @@ type t = {
center_a : Coordinate.t;
center : Coordinate.t;
norm_sq : float;
norm : float;
norm_coef: float;
coef : float;
norm_fun : int array -> int array -> float;
i : int;
@ -31,28 +31,39 @@ let create_array ?cutoff p_a p_b =
let norm_sq =
Coordinate.dot center_ab center_ab
in
let norm_coef_scale_a =
Contracted_shell.norm_coef_scale p_a
and norm_coef_scale_b =
Contracted_shell.norm_coef_scale p_b
in
let norm_fun a b =
let k1, k2 =
Zkey.of_int_array a ~kind:Kind_3,
Zkey.of_int_array b ~kind:Kind_3
in
let v1 =
Zmap.find norm_coef_scale_a k1
and v2 =
Zmap.find norm_coef_scale_b k2
in v1 *. v2
in
Array.init (Contracted_shell.size p_a) (fun i ->
let p_a_expo_center = Coordinate.(
Contracted_shell.expo p_a i |. Contracted_shell.center p_a )
in
let f1 =
let norm_coef_a =
Contracted_shell.norm_coef p_a i
in
Array.init (Contracted_shell.size p_b) (fun j ->
try
let f2 =
let norm_coef_b =
Contracted_shell.norm_coef p_b j
in
let norm_fun a b =
f1 a *. f2 b
let norm_coef =
norm_coef_a *. norm_coef_b
in
let norm =
norm_fun
[| Angular_momentum.to_int @@ Contracted_shell.totAngMom p_a ; 0 ; 0 |]
[| Angular_momentum.to_int @@ Contracted_shell.totAngMom p_b ; 0 ; 0 |]
in
if (norm < cutoff) then
if (norm_coef < cutoff) then
raise Null_contribution;
let p_b_expo_center = Coordinate.(
Contracted_shell.expo p_b j |. Contracted_shell.center p_b )
@ -70,19 +81,15 @@ let create_array ?cutoff p_a p_b =
let g =
(pi *. expo_inv)**(1.5) *. exp(-. argexpo)
in
let norm_inv = 1./.norm in
let norm_fun a b =
norm_inv *. norm_fun a b
in
let coef =
norm *. Contracted_shell.(coef p_a i *. coef p_b j) *. g
norm_coef *. Contracted_shell.(coef p_a i *. coef p_b j) *. g
in
if (abs_float coef < cutoff) then
raise Null_contribution;
let center_a =
Coordinate.(center |- Contracted_shell.center p_a)
in
Some { i ; j ; shell_a=p_a ; shell_b=p_b ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq }
Some { i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; norm_fun ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq }
with
| Null_contribution -> None
)

215
Basis/TwoElectronRR.ml
View File

@ -1,6 +1,7 @@
open Util
let cutoff2 = cutoff *. cutoff
let debug = false
exception NullQuartet
@ -22,78 +23,97 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let maxm = totAngMom_a + totAngMom_b + totAngMom_c + totAngMom_d in
let empty = Array.make (maxm+1) 0.
in
if debug then begin
Printf.printf "\n---- %d %d %d %d ----\n" totAngMom_a totAngMom_b totAngMom_c totAngMom_d;
Printf.printf "%d %d %d\n" angMom_a.(0) angMom_a.(1) angMom_a.(2) ;
Printf.printf "%d %d %d\n" angMom_b.(0) angMom_b.(1) angMom_b.(2) ;
Printf.printf "%d %d %d\n" angMom_c.(0) angMom_c.(1) angMom_c.(2) ;
Printf.printf "%d %d %d\n" angMom_d.(0) angMom_d.(1) angMom_d.(2) ;
Printf.printf "%f %f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n" expo_b expo_d
end
expo_inv_p expo_inv_q
(Coordinate.coord center_ab 0) (Coordinate.coord center_ab 1) (Coordinate.coord center_ab 2)
(Coordinate.coord center_cd 0) (Coordinate.coord center_cd 1) (Coordinate.coord center_cd 2)
(Coordinate.coord center_pq 0) (Coordinate.coord center_pq 1) (Coordinate.coord center_pq 2);
(** Vertical recurrence relations *)
let rec vrr0 angMom_a = function
let rec vrr0 angMom_a totAngMom_a =
if debug then
Printf.printf "vrr0: %d : %d %d %d\n" totAngMom_a angMom_a.(0) angMom_a.(1) angMom_a.(2);
match totAngMom_a with
| 0 -> zero_m_array
| totAngMom_a ->
| _ ->
let maxsze = maxm+1 in
let key = Zkey.of_int_tuple (Zkey.Three
(angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1) )
in
let (found, result) =
try (true, Zmap.find map key) with
| Not_found -> (false,
let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and amm = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and xyz =
match angMom_a with
| [|0;0;_|] -> 2
| [|0;_;_|] -> 1
| _ -> 0
in
am.(xyz) <- am.(xyz) - 1;
amm.(xyz) <- amm.(xyz) - 2;
if am.(xyz) < 0 then empty else
let v1 =
vrr0 am (totAngMom_a-1)
in
let f1 = expo_inv_p *. (Coordinate.coord center_pq xyz)
and f2 = expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)
in
if amm.(xyz) < 0 then
Array.init (maxm+1) (fun m ->
if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
else
let f3 = (float_of_int am.(xyz)) *. expo_inv_p *. 0.5 in
let v3 =
vrr0 amm (totAngMom_a-2)
in
Array.init (maxm+1) (fun m ->
(if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
+. f3 *. (v3.(m) +. if m = maxm then 0. else
expo_inv_p *. v3.(m+1))
)
)
in
if not found then
Zmap.add map key result;
try Zmap.find map key with
| Not_found ->
let result =
let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and amm = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
in
let xyz =
match angMom_a with
| [|_;0;0|] -> 0
| [|_;_;0|] -> 1
| _ -> 2
in
am.(xyz) <- am.(xyz) - 1;
amm.(xyz) <- amm.(xyz) - 2;
if am.(xyz) < 0 then empty else
let v1 =
vrr0 am (totAngMom_a-1)
in
let f1 = expo_inv_p *. (Coordinate.coord center_pq xyz)
and f2 = expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)
in
if amm.(xyz) < 0 then
Array.init (maxsze) (fun m ->
if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
else
let f3 = (float_of_int am.(xyz)) *. expo_inv_p *. 0.5 in
let v3 =
vrr0 amm (totAngMom_a-2)
in
Array.init (maxsze) (fun m ->
(if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) )
+. f3 *. (v3.(m) +. if m = maxm then 0. else
expo_inv_p *. v3.(m+1))
)
in Zmap.add map key result;
result
and vrr angMom_a angMom_c totAngMom_a totAngMom_c =
if debug then
Printf.printf "vrr : %d %d : %d %d %d %d %d %d\n" totAngMom_a totAngMom_c angMom_a.(0) angMom_a.(1) angMom_a.(2) angMom_c.(0) angMom_c.(1) angMom_c.(2);
match (totAngMom_a, totAngMom_c) with
| (0,0) -> zero_m_array
| (_,0) -> vrr0 angMom_a totAngMom_a
| (_,_) ->
let maxsze = maxm+1 in
let key = Zkey.of_int_tuple (Zkey.Six
((angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1),
(angMom_c.(0)+1, angMom_c.(1)+1, angMom_c.(2)+1)) )
in
let (found, result) =
try (true, Zmap.find map key) with
| Not_found -> (false,
try Zmap.find map key with
| Not_found ->
let result =
let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and cm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |]
and cmm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |]
and xyz =
match angMom_c with
| [|0;0;_|] -> 2
| [|0;_;_|] -> 1
| _ -> 0
| [|_;0;0|] -> 0
| [|_;_;0|] -> 1
| _ -> 2
in
am.(xyz) <- am.(xyz) - 1;
cm.(xyz) <- cm.(xyz) - 1;
@ -110,7 +130,7 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let v1 =
vrr angMom_a cm totAngMom_a (totAngMom_c-1)
in
Array.init (maxm+1) (fun m ->
Array.init (maxsze) (fun m ->
f1 *. v1.(m) -. (if m = maxm then 0. else f2 *. v1.(m+1)) )
in
let result =
@ -118,11 +138,11 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let f3 =
(float_of_int cm.(xyz)) *. expo_inv_q *. 0.5
in
if (abs_float f3 < cutoff) && (abs_float (f3 *. abs_float expo_inv_q) < cutoff) then result else
if (abs_float f3 < cutoff) && (abs_float (f3 *. expo_inv_q) < cutoff) then result else
let v3 =
vrr angMom_a cmm totAngMom_a (totAngMom_c-2)
in
Array.init (maxm+1) (fun m -> result.(m) +.
Array.init (maxsze) (fun m -> result.(m) +.
f3 *. (v3.(m) +. (if m=maxm then 0. else expo_inv_q *. v3.(m+1)) ))
in
let result =
@ -134,66 +154,72 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
let v5 =
vrr am cm (totAngMom_a-1) (totAngMom_c-1)
in
Array.init (maxm+1) (fun m ->
Array.init (maxsze) (fun m ->
result.(m) -. (if m = maxm then 0. else f5 *. v5.(m+1)))
in
result
)
in
if not found then
Zmap.add map key result;
in Zmap.add map key result;
result
(** Horizontal recurrence relations *)
and hrr0 angMom_a angMom_b angMom_c
totAngMom_a totAngMom_b totAngMom_c =
if debug then
Printf.printf "hrr0: %d %d %d : %d %d %d %d %d %d %d %d %d\n" totAngMom_a totAngMom_b totAngMom_c angMom_a.(0) angMom_a.(1) angMom_a.(2) angMom_b.(0) angMom_b.(1) angMom_b.(2) angMom_c.(0) angMom_c.(1) angMom_c.(2);
match totAngMom_b with
| 0 -> (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0)
| 1 -> let xyz = if angMom_b.(0) = 1 then 0 else if angMom_b.(1) = 1 then 1 else 2 in
let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in
ap.(xyz) <- ap.(xyz) + 1;
let v1 =
vrr ap angMom_c (totAngMom_a+1) totAngMom_c
in
let f2 =
(Coordinate.coord center_ab xyz)
in
if (abs_float f2 < cutoff) then v1.(0) else
let v2 =
vrr angMom_a angMom_c totAngMom_a totAngMom_c
in
v1.(0) +. f2 *. v2.(0)
| 1 ->
let xyz = if angMom_b.(0) = 1 then 0 else if angMom_b.(1) = 1 then 1 else 2 in
let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in
ap.(xyz) <- ap.(xyz) + 1;
let v1 =
vrr ap angMom_c (totAngMom_a+1) totAngMom_c
in
let f2 =
(Coordinate.coord center_ab xyz)
in
if (abs_float f2 < cutoff) then v1.(0) else
let v2 =
vrr angMom_a angMom_c totAngMom_a totAngMom_c
in
v1.(0) +. f2 *. v2.(0)
| _ ->
let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |]
and xyz =
match angMom_b with
| [|0;0;_|] -> 2
| [|0;_;_|] -> 1
| _ -> 0
in
ap.(xyz) <- ap.(xyz) + 1;
bm.(xyz) <- bm.(xyz) - 1;
if (bm.(xyz) < 0) then 0. else
let h1 =
hrr0 ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c
let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |]
and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |]
and xyz =
match angMom_b with
| [|_;0;0|] -> 0
| [|_;_;0|] -> 1
| _ -> 2
in
let f2 =
(Coordinate.coord center_ab xyz)
in
if (abs_float f2 < cutoff) then h1 else
let h2 =
hrr0 angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c
ap.(xyz) <- ap.(xyz) + 1;
bm.(xyz) <- bm.(xyz) - 1;
if (bm.(xyz) < 0) then 0. else
let h1 =
hrr0 ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c
in
h1 +. f2 *. h2
let f2 =
(Coordinate.coord center_ab xyz)
in
if (abs_float f2 < cutoff) then h1 else
let h2 =
hrr0 angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c
in
h1 +. f2 *. h2
and hrr angMom_a angMom_b angMom_c angMom_d
totAngMom_a totAngMom_b totAngMom_c totAngMom_d =
if debug then
Printf.printf "hrr : %d %d %d %d : %d %d %d %d %d %d %d %d %d %d %d %d\n" totAngMom_a totAngMom_b totAngMom_c totAngMom_d angMom_a.(0) angMom_a.(1) angMom_a.(2) angMom_b.(0) angMom_b.(1) angMom_b.(2) angMom_c.(0) angMom_c.(1) angMom_c.(2) angMom_d.(0) angMom_d.(1) angMom_d.(2);
match (totAngMom_b, totAngMom_d) with
| (0,0) -> (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0)
| (_,0) -> hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c
@ -202,9 +228,9 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
and dm = [| angMom_d.(0) ; angMom_d.(1) ; angMom_d.(2) |]
and xyz =
match angMom_d with
| [|0;0;_|] -> 2
| [|0;_;_|] -> 1
| _ -> 0
| [|_;0;0|] -> 0
| [|_;_;0|] -> 1
| _ -> 2
in
cp.(xyz) <- cp.(xyz) + 1;
dm.(xyz) <- dm.(xyz) - 1;
@ -217,9 +243,9 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d)
hrr angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1)
in
h1 +. f2 *. h2
in
hrr angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b
totAngMom_c totAngMom_d
hrr angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d
@ -288,7 +314,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
let d = shell_q.(cd).Shell_pair.j in
let map = Zmap.create (Array.length class_indices) in
(* Compute the integral class from the primitive shell quartet *)
Array.iteri (fun i key ->
class_indices
|> Array.iteri (fun i key ->
let a = Zkey.to_int_array Zkey.Kind_12 key in
let (angMomA,angMomB,angMomC,angMomD) =
( [| a.(0) ; a.(1) ; a.(2) |],
@ -343,7 +370,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
in
contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral
with NullQuartet -> ()
) class_indices
)
end
with NullQuartet -> ()
done

2
Makefile
View File

@ -42,4 +42,4 @@ clean:
rm -rf _build $(ALL_EXE) $(ALL_TESTS) *.native *.byte
debug: run_integrals.native
time ./run_integrals -c h2o.xyz -b ~/quantum_package/data/basis/cc-pvtz -o /dev/shm/out ; sleep 2 ; diff /dev/shm/out.eri REF | head -50
./debug.sh