diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml index 990e586..47c778a 100644 --- a/Basis/Contracted_shell.ml +++ b/Basis/Contracted_shell.ml @@ -7,7 +7,7 @@ type t = { totAngMom : Angular_momentum.t; size : int; norm_coef : float array; - norm_coef_scale : float Zmap.t; + norm_coef_scale : float array; indice : int; powers : Zkey.t array; } @@ -82,10 +82,11 @@ let create ~indice ~expo ~coef ~center ~totAngMom = 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 + let norm_coef_scale = + Array.map (fun a -> + (norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0) + ) powers 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 ; norm_coef_scale ; powers } diff --git a/Basis/ERI.ml b/Basis/ERI.ml index c69d9a4..27776fa 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -93,6 +93,7 @@ let to_file ~filename basis = 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 @@ -102,6 +103,7 @@ let to_file ~filename basis = let shell_q = shell_pairs.(k).(l) in + let swap = Array.length shell_q < Array.length shell_p in @@ -144,6 +146,9 @@ let to_file ~filename basis = in if (abs_float value > cutoff) then (inn := !inn + 1; + (* + assert ((i_c, k_c, j_c, l_c) <> (-1,0,0,0)); + *) Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" i_c k_c j_c l_c value ) diff --git a/Basis/Shell_pair.ml b/Basis/Shell_pair.ml index 96ff176..f39f171 100644 --- a/Basis/Shell_pair.ml +++ b/Basis/Shell_pair.ml @@ -9,7 +9,7 @@ type t = { norm_sq : float; norm_coef: float; coef : float; - norm_fun : int array -> int array -> float; + norm_coef_scale : float array; i : int; j : int; shell_a : Contracted_shell.t; @@ -36,16 +36,12 @@ let create_array ?cutoff p_a p_b = 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 + let norm_coef_scale = + Array.map (fun v1 -> + Array.map (fun v2 -> v1 *. v2) norm_coef_scale_b + ) norm_coef_scale_a + |> Array.to_list + |> Array.concat in Array.init (Contracted_shell.size p_a) (fun i -> let p_a_expo_center = Coordinate.( @@ -89,7 +85,7 @@ let create_array ?cutoff p_a p_b = let center_a = Coordinate.(center |- Contracted_shell.center p_a) in - 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 } + Some { i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; norm_coef_scale } with | Null_contribution -> None ) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index eb87d4c..35b6b29 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -280,6 +280,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q 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 + let norm_coef_scale_p = shell_p.(ab).Shell_pair.norm_coef_scale in for cd=0 to (Array.length shell_q - 1) do let coef_prod = cab *. shell_q.(cd).Shell_pair.coef @@ -313,6 +314,14 @@ 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 + let norm_coef_scale_q = shell_q.(cd).Shell_pair.norm_coef_scale in + let norm_coef_scale = + Array.map (fun v1 -> + Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q + ) norm_coef_scale_p + |> Array.to_list + |> Array.concat + in (* Compute the integral class from the primitive shell quartet *) class_indices |> Array.iteri (fun i key -> @@ -355,9 +364,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q *) - let norm = - shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD - in + let norm = norm_coef_scale.(i) in let integral = chop norm (fun () -> hvrr_two_e (angMomA, angMomB, angMomC, angMomD) (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 8ca3bdf..d813183 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -324,13 +324,14 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q accu +. coef_prod *. zero_m_array.(0) with NullQuartet -> accu ) 0. shell_q - ) 0. shell_p + ) 0. shell_p | _ -> Array.iter (fun shell_ab -> + let norm_coef_scale_p = shell_ab.Shell_pair.norm_coef_scale in let b = shell_ab.Shell_pair.j in let common = - Array.mapi (fun idx shell_cd -> + Array.map (fun shell_cd -> let coef_prod = shell_ab.Shell_pair.coef *. shell_cd.Shell_pair.coef in @@ -352,30 +353,36 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (zero_m_array, shell_cd.Shell_pair.expo_inv, Contracted_shell.expo shell_d d, shell_cd.Shell_pair.center_ab, - center_pq,coef_prod,idx) + center_pq,coef_prod) ) shell_q |> Array.to_list |> List.filter (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod,idx) -> abs_float coef_prod >= 1.e-4 *. cutoff) + center_pq,coef_prod) -> abs_float coef_prod >= 1.e-4 *. cutoff) |> Array.of_list in let zero_m_array = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod,idx) -> zero_m_array) common + center_pq,coef_prod) -> zero_m_array) common and expo_inv = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod,idx) -> expo_inv ) common + center_pq,coef_prod) -> expo_inv ) common and d = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod,idx) -> d) common + center_pq,coef_prod) -> d) common and center_cd = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod,idx) -> center_cd) common + center_pq,coef_prod) -> center_cd) common and center_pq = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod,idx) -> center_pq) common + center_pq,coef_prod) -> center_pq) common and coef_prod = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod,idx) -> coef_prod) common - and idx = Array.map (fun (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod,idx) -> idx) common + center_pq,coef_prod) -> coef_prod) common in (* Compute the integral class from the primitive shell quartet *) let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in + let norm = + let norm_coef_scale_q = shell_q.(0).Shell_pair.norm_coef_scale in + Array.map (fun v1 -> + Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q + ) norm_coef_scale_p + |> Array.to_list + |> Array.concat + in Array.iteri (fun i key -> let a = Zkey.to_int_array Zkey.Kind_12 key in let (angMomA,angMomB,angMomC,angMomD) = @@ -384,11 +391,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q [| a.(6) ; a.(7) ; a.(8) |], [| a.(9) ; a.(10) ; a.(11) |] ) in - let norm = - Array.map (fun shell_cd -> - shell_ab.Shell_pair.norm_fun angMomA angMomB *. shell_cd.Shell_pair.norm_fun angMomC angMomD - ) shell_q - in let integral = hvrr_two_e_vector (angMomA, angMomB, angMomC, angMomD) (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, @@ -398,7 +400,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (shell_ab.Shell_pair.expo_inv, expo_inv) (shell_ab.Shell_pair.center_ab, center_cd, center_pq) coef_prod map - |> Array.mapi (fun i x -> x *. norm.(idx.(i)) ) + |> Array.map (fun x -> x *. norm.(i) ) in let x = Array.fold_left (+.) 0. integral in contracted_class.(i) <- contracted_class.(i) +. x