diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index ddb0891..3a4f82d 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -129,17 +129,17 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let angMom_ax, angMom_ay, angMom_az = angMom_a and angMom_cx, angMom_cy, angMom_cz = angMom_c in match angMom_c with - | (_,0,0) -> (* 321_984 *) + | (_,0,0) -> (angMom_ax-1, angMom_ay, angMom_az), (angMom_cx-1, angMom_cy, angMom_cz), (angMom_cx-2, angMom_cy, angMom_cz), angMom_ax,angMom_cx-1, 0 - | (_,_,0) -> (* 612_002 *) + | (_,_,0) -> (angMom_ax, angMom_ay-1, angMom_az), (angMom_cx, angMom_cy-1, angMom_cz), (angMom_cx, angMom_cy-2, angMom_cz), angMom_ay,angMom_cy-1, 1 - | _ -> (* 1_067_324 *) + | _ -> (angMom_ax, angMom_ay, angMom_az-1), (angMom_cx, angMom_cy, angMom_cz-1), (angMom_cx, angMom_cy, angMom_cz-2), @@ -152,14 +152,13 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let f2 = expo_inv_q *. (Coordinate.coord center_pq xyz) in - let result = empty in let result = - if ( (abs_float f1 < cutoff) && (abs_float f2 < cutoff) ) then result else + if ( (abs_float f1 < cutoff) && (abs_float f2 < cutoff) ) then empty else let v1 = vrr angMom_a cm totAngMom_a (totAngMom_c-1) in - Array.init maxsze (fun m -> result.(m) +. - f1 *. v1.(m) -. (if m = maxm then 0. else f2 *. v1.(m+1)) ) + Array.init maxsze (fun m -> f1 *. v1.(m) -. + (if m = maxm then 0. else f2 *. v1.(m+1)) ) in let result = if cmxyz < 1 then result else diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 29453f5..346615f 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -497,6 +497,61 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.make (Array.length class_indices) 0.; in + (** Screening on the product of coefficients *) + let coef_max_p = + Array.fold_left (fun accu x -> + if (abs_float x) > accu then (abs_float x) else accu) + 0. shell_p.ContractedShellPair.coef + and coef_max_q = + Array.fold_left (fun accu x -> + if (abs_float x) > accu then (abs_float x) else accu) + 0. shell_q.ContractedShellPair.coef + in + + let rec build_list cutoff vec accu = function + | -1 -> Array.of_list accu + | k -> build_list cutoff vec ( + if (abs_float vec.(k) > cutoff) then (k::accu) + else accu ) (k-1) + in + let p_list = + let vec = shell_p.ContractedShellPair.coef in + build_list (cutoff /. coef_max_q) vec [] (Array.length vec - 1) + and q_list = + let vec = shell_q.ContractedShellPair.coef in + build_list (cutoff /. coef_max_p) vec [] (Array.length vec - 1) + in + +(* + if (Array.length p_list) < (Array.length shell_p.ContractedShellPair.coef) then + begin + Printf.printf "Reduced p from %d to %d:\n" (Array.length + shell_p.ContractedShellPair.coef) (Array.length p_list); + Array.iter (fun k -> Printf.printf "%d " k) p_list; print_newline (); + Printf.printf "\n%!" + end; + + if (Array.length q_list) < (Array.length shell_q.ContractedShellPair.coef) then + begin + Printf.printf "Reduced q from %d to %d:\n" (Array.length + shell_q.ContractedShellPair.coef) (Array.length q_list); + Array.iter (fun k -> Printf.printf "%d " k) q_list; print_newline (); + Printf.printf "\n%!" + end; + *) + + let np, nq = + Array.length p_list, + Array.length q_list + in + let filter_p vec = Array.init np (fun k -> vec.(p_list.(k))) + and filter_q vec = Array.init nq (fun k -> vec.(q_list.(k))) + in + let sp = filter_p sp + and sq = filter_q sq + in + + (* Compute all integrals in the shell for each pair of significant shell pairs *) begin @@ -504,9 +559,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q totAngMom shell_c, totAngMom shell_d) with | Angular_momentum.(S,S,S,S) -> contracted_class.(0) <- - let np, nq = - Array.length sp, Array.length sq - in + begin + try let expo_inv_p = Vec.init np (fun ab -> sp.(ab-1).ShellPair.expo_inv) and expo_inv_q = @@ -516,13 +570,12 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let coef = let result = Mat.make0 nq np in Lacaml.D.ger - (Vec.of_array shell_q.ContractedShellPair.coef) - (Vec.of_array shell_p.ContractedShellPair.coef) + (Vec.of_array @@ filter_q shell_q.ContractedShellPair.coef) + (Vec.of_array @@ filter_p shell_p.ContractedShellPair.coef) result; result in let zm_array = Mat.init_cols np nq (fun i j -> - (** Screening on the product of coefficients *) try if (abs_float coef.{j,i} ) < 1.e-3*.cutoff then raise NullQuartet; @@ -544,15 +597,23 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q with NullQuartet -> 0. ) in Mat.gemm_trace zm_array coef + with (Invalid_argument _) -> 0. + end | _ -> + + let coef = + let cp = filter_p shell_p.ContractedShellPair.coef + and cq = filter_q shell_q.ContractedShellPair.coef + in + Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) ) + in + let expo_inv_p = - Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp + Array.map (fun shell_ab -> shell_ab.ShellPair.expo_inv) sp and expo_inv_q = - Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq - in - let np, nq = - Array.length expo_inv_p, Array.length expo_inv_q + Array.map (fun shell_cd -> shell_cd.ShellPair.expo_inv) sq in + let expo_b = Array.map (fun shell_ab -> Contracted_shell.expo shell_b shell_ab.ShellPair.j) sp and expo_d = @@ -560,14 +621,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in - let coef = - Array.init np (fun l -> - Array.init nq (fun k -> - shell_q.ContractedShellPair.coef.(k) *. - shell_p.ContractedShellPair.coef.(l) - ) ) - in - let center_pq = Array.init 3 (fun xyz -> Array.init np (fun ab -> @@ -610,12 +663,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q zero_m ~maxm ~expo_pq_inv ~norm_pq_sq ) sq - (* - |> Array.to_list - |> List.filter (fun (zero_m_array, d, - center_pq,coef_prod) -> abs_float coef_prod >= 1.e-4 *. cutoff) - |> Array.of_list - *) in (* Transpose result *) for m=0 to maxm do