From f3d03cd4246d660e4cb68074ed2f986f537a75b3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 1 Feb 2018 16:09:04 +0100 Subject: [PATCH] Accelerated norm function --- Basis/Contracted_shell.ml | 18 +++- Basis/Shell_pair.ml | 41 +++++--- Basis/TwoElectronRR.ml | 215 +++++++++++++++++++++----------------- Makefile | 2 +- 4 files changed, 161 insertions(+), 115 deletions(-) diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml index 798fabd..990e586 100644 --- a/Basis/Contracted_shell.ml +++ b/Basis/Contracted_shell.ml @@ -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 } diff --git a/Basis/Shell_pair.ml b/Basis/Shell_pair.ml index d7fb4a7..96ff176 100644 --- a/Basis/Shell_pair.ml +++ b/Basis/Shell_pair.ml @@ -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 ) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 773512d..eb87d4c 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -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 diff --git a/Makefile b/Makefile index 71911b6..92b2ee1 100644 --- a/Makefile +++ b/Makefile @@ -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