From e0c0ec13535bc814a918c2f6f89f64c2016ba64b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 23 Jan 2018 19:26:28 +0100 Subject: [PATCH] Schwartz --- Basis/ERI.ml | 103 ++++++++++++++++++++++++++++------------- Basis/Shell_pair.ml | 4 +- Basis/TwoElectronRR.ml | 68 +++++++++++++++++---------- 3 files changed, 117 insertions(+), 58 deletions(-) diff --git a/Basis/ERI.ml b/Basis/ERI.ml index bb6f559..93a1cb6 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -25,8 +25,13 @@ let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t = TwoElectronRR.contracted_class ~zero_m shell_a shell_b shell_c shell_d +(** Compute all the integrals of a contracted class *) +let contracted_class_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = + TwoElectronRR.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q + type n_cls = { n : int ; cls : Z.t array } +exception NullIntegral (** Write all integrals to a file with the convention *) let to_file ~filename basis = @@ -38,44 +43,76 @@ let to_file ~filename basis = in let oc = open_out filename in + + (* Pre-compute all shell pairs *) + let shell_pairs = + Array.mapi (fun i shell_a -> Array.map (fun shell_b -> + Shell_pair.create_array shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis + in + + (* Pre-compute diagonal integrals for Schwartz *) + let schwartz = + Array.map (fun pair_array -> Array.map (fun pair -> + let cls = + contracted_class_shell_pairs pair pair + in + (cls, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) + ) pair_array ) shell_pairs + in + for i=0 to (Array.length basis) - 1 do print_int basis.(i).Contracted_shell.indice ; print_newline (); for j=0 to i do - for k=0 to i do - for l=0 to k do - (* Compute all the integrals of the class *) - let cls = - contracted_class basis.(i) basis.(j) basis.(k) basis.(l) - in + let schwartz_p, schwartz_p_max = schwartz.(i).(j) in + try + if (schwartz_p_max < cutoff) then raise NullIntegral; + let + shell_p = shell_pairs.(i).(j) + in + for k=0 to i do + for l=0 to k do + let schwartz_q, schwartz_q_max = schwartz.(k).(l) in + try + if schwartz_p_max *. schwartz_q_max < cutoff *. cutoff then + raise NullIntegral; + let + shell_q = shell_pairs.(k).(l) + in + (* Compute all the integrals of the class *) + let cls = + contracted_class_shell_pairs ~schwartz_p ~schwartz_q shell_p shell_q + in - (* Write the data in the output file *) - Array.iteri (fun i_c powers_i -> - let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in - let xi = to_int_tuple powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in - let xj = to_int_tuple powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = basis.(k).Contracted_shell.indice + k_c + 1 in - let xk = to_int_tuple powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = basis.(l).Contracted_shell.indice + l_c + 1 in - let xl = to_int_tuple powers_l in - let key = - Zkey.of_int_tuple (Zkey.Twelve (xi,xj,xk,xl)) - in - let value = - Zmap.find cls key - in - if (abs_float value > cutoff) then - Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" - i_c k_c j_c l_c value - ) basis.(l).Contracted_shell.powers - ) basis.(k).Contracted_shell.powers - ) basis.(j).Contracted_shell.powers - ) basis.(i).Contracted_shell.powers; + (* Write the data in the output file *) + Array.iteri (fun i_c powers_i -> + let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in + let xi = to_int_tuple powers_i in + Array.iteri (fun j_c powers_j -> + let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in + let xj = to_int_tuple powers_j in + Array.iteri (fun k_c powers_k -> + let k_c = basis.(k).Contracted_shell.indice + k_c + 1 in + let xk = to_int_tuple powers_k in + Array.iteri (fun l_c powers_l -> + let l_c = basis.(l).Contracted_shell.indice + l_c + 1 in + let xl = to_int_tuple powers_l in + let key = + Zkey.of_int_tuple (Zkey.Twelve (xi,xj,xk,xl)) + in + let value = + Zmap.find cls key + in + if (abs_float value > cutoff) then + Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" + i_c k_c j_c l_c value + ) basis.(l).Contracted_shell.powers + ) basis.(k).Contracted_shell.powers + ) basis.(j).Contracted_shell.powers + ) basis.(i).Contracted_shell.powers; + with NullIntegral -> print_endline "Schwartz"; () + done; done; - done; + with NullIntegral -> print_endline "Big Schwartz"; () done; done; close_out oc diff --git a/Basis/Shell_pair.ml b/Basis/Shell_pair.ml index 240604d..d7fb4a7 100644 --- a/Basis/Shell_pair.ml +++ b/Basis/Shell_pair.ml @@ -12,6 +12,8 @@ type t = { norm_fun : int array -> int array -> float; i : int; j : int; + shell_a : Contracted_shell.t; + shell_b : Contracted_shell.t; } exception Null_contribution @@ -80,7 +82,7 @@ let create_array ?cutoff p_a p_b = let center_a = Coordinate.(center |- Contracted_shell.center p_a) in - Some { i ; j ; 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_fun ; norm ; 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 32de469..1338ba4 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -168,13 +168,14 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d) +let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = -(** Computes all the two-electron integrals of the contracted shell quartet *) -let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = - - let shell_p = Shell_pair.create_array shell_a shell_b - and shell_q = Shell_pair.create_array shell_c shell_d - and maxm = + let shell_a = shell_p.(0).Shell_pair.shell_a + and shell_b = shell_p.(0).Shell_pair.shell_b + and shell_c = shell_q.(0).Shell_pair.shell_a + and shell_d = shell_q.(0).Shell_pair.shell_b + in + let maxm = let open Angular_momentum in (to_int @@ Contracted_shell.totAngMom shell_a) + (to_int @@ Contracted_shell.totAngMom shell_b) + (to_int @@ Contracted_shell.totAngMom shell_c) + (to_int @@ Contracted_shell.totAngMom shell_d) @@ -248,31 +249,44 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in (* Compute the integral class from the primitive shell quartet *) Array.iteri (fun i key -> + let a = Zkey.to_int_array Zkey.Kind_12 key in let (angMomA,angMomB,angMomC,angMomD) = - let a = Zkey.to_int_array Zkey.Kind_12 key in ( [| a.(0) ; a.(1) ; a.(2) |], [| a.(3) ; a.(4) ; a.(5) |], [| a.(6) ; a.(7) ; a.(8) |], [| a.(9) ; a.(10) ; a.(11) |] ) in try - (* (* Schwartz screening *) - let schwartz_p = - Zmap.find overlaps_p @@ Zkey.of_int_array ~kind:Zkey.Kind_6 - [| angMomA.(0) ; angMomA.(1) ; angMomA.(2) ; - angMomB.(0) ; angMomB.(1) ; angMomB.(2) |] - |> abs_float + (* + let schwartz_p = + let key = + Zkey.of_int_array Zkey.Kind_12 + [| a.(0) ; a.(1) ; a.(2) ; + a.(3) ; a.(4) ; a.(5) ; + a.(0) ; a.(1) ; a.(2) ; + a.(3) ; a.(4) ; a.(5) |] + in + match schwartz_p with + | None -> 1. + | Some schwartz_p -> Zmap.find schwartz_p key in - let schwartz_q = - Zmap.find overlaps_q @@ Zkey.of_int_array ~kind:Zkey.Kind_6 - [| angMomC.(0) ; angMomC.(1) ; angMomC.(2) ; - angMomD.(0) ; angMomD.(1) ; angMomD.(2) |] - |> abs_float + if schwartz_p < cutoff then raise NullQuartet; + let schwartz_q = + let key = + Zkey.of_int_array Zkey.Kind_12 + [| a.(6) ; a.(7) ; a.(8) ; + a.(9) ; a.(10) ; a.(11) ; + a.(6) ; a.(7) ; a.(8) ; + a.(9) ; a.(10) ; a.(11) |] + in + match schwartz_q with + | None -> 1. + | Some schwartz_q -> Zmap.find schwartz_q key in - if schwartz_p*.schwartz_q = 0. then - () ; (*raise NullQuartet; *) - *) + if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet; + *) + let norm = shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD @@ -288,9 +302,6 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = map ) in contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral - (* -;if (schwartz_p*.schwartz_q < cutoff2) then Printf.printf "%e %e\n" (schwartz_p*.schwartz_q) integral; - *) with NullQuartet -> () ) class_indices with NullQuartet -> () @@ -304,3 +315,12 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = +(** Computes all the two-electron integrals of the contracted shell quartet *) +let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = + + let shell_p = Shell_pair.create_array shell_a shell_b + and shell_q = Shell_pair.create_array shell_c shell_d + in + contracted_class_shell_pairs ~zero_m shell_p shell_q + +