diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index 68590c3..99f48ab 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -16,7 +16,8 @@ type t = center_ab : Coordinate.t; (* A-B *) norm_sq : float; (* |A-B|^2 *) norm_coef_scale : float array; (* norm_coef.(i) / norm_coef.(0) *) - totAngMomInt : int (* Total angular Momentum *) + totAngMomInt : int; (* Total angular Momentum *) + monocentric : bool; } @@ -119,6 +120,7 @@ let create ?cutoff p_a p_b = shell_pairs ; center_ab=shell_pairs.(0).center_ab; norm_coef_scale ; norm_sq=shell_pairs.(0).norm_sq; totAngMomInt = shell_pairs.(0).ShellPair.totAngMomInt; + monocentric = shell_pairs.(0).ShellPair.monocentric; } diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index e620999..da5f37f 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -86,11 +86,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) done; result.(maxm) <- -. f2 *. v1.(maxm) - (* - Array.init maxsze (fun m -> - if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) ) - *) - end + end else begin let v3 = @@ -107,14 +103,6 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) result.(maxm) <- f3 *. v3.(maxm) end; result - (* - 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_1d key result; result @@ -164,47 +152,55 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (angMom_cx, angMom_cy, angMom_cz-2), angMom_az,angMom_cz-1, 2 in - if cmxyz < 0 then empty else + if cmxyz < 0 then empty + else let f1 = -. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz) - in - let f2 = + and f2 = expo_inv_q *. (Coordinate.coord center_pq xyz) in - let result = - if ( (abs_float f1 < cutoff) && (abs_float f2 < cutoff) ) then empty else + let result = Array.make maxsze 0. in + if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then + begin let v1 = vrr angMom_a cm totAngMom_a (totAngMom_c-1) in - 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 + for m=0 to maxm-1 do + result.(m) <- f1 *. v1.(m) -. f2 *. v1.(m+1) ; + done; + result.(maxm) <- f1 *. v1.(maxm) ; + end; + if cmxyz > 0 then + begin let f3 = (float_of_int cmxyz) *. expo_inv_q *. 0.5 in - if (abs_float f3 < cutoff) && (abs_float (f3 *. expo_inv_q) < cutoff) then result else - ( + if (abs_float f3 > cutoff) || + (abs_float (f3 *. expo_inv_q) > cutoff) then + begin let v3 = vrr angMom_a cmm totAngMom_a (totAngMom_c-2) in - 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 = - if (axyz < 1) || (cmxyz < 0) then result else + for m=0 to maxm-1 do + result.(m) <- result.(m) +. + f3 *. (v3.(m) +. expo_inv_q *. v3.(m+1)) + done; + result.(maxm) <- result.(maxm) +. f3 *. v3.(maxm) + end + end; + if (axyz > 0) && (cmxyz >= 0) then + begin let f5 = (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 in - if (abs_float f5 < cutoff) then result else + if (abs_float f5 > cutoff) then let v5 = vrr am cm (totAngMom_a-1) (totAngMom_c-1) in - Array.init (maxsze) (fun m -> - result.(m) -. (if m = maxm then 0. else f5 *. v5.(m+1))) - in + for m=0 to maxm-1 do + result.(m) <- result.(m) -. f5 *. v5.(m+1) + done + end; result in Zmap.add map_2d key result; result @@ -352,6 +348,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.make (Array.length class_indices) 0.; in + let monocentric = + shell_p.ContractedShellPair.monocentric && + shell_q.ContractedShellPair.monocentric + in + (* Compute all integrals in the shell for each pair of significant shell pairs *) let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in @@ -365,7 +366,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in (** Screening on the product of coefficients *) try - if (abs_float coef_prod) < 1.e-4*.cutoff then + if (abs_float coef_prod) < 1.e-3*.cutoff then raise NullQuartet; let expo_pq_inv = @@ -397,18 +398,19 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let map_1d' = Zmap.create (4*maxm) in let map_2d' = Zmap.create (Array.length class_indices) 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.to_list norm_coef_scale_p + |> List.map (fun v1 -> + Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) |> Array.concat in - (* - let monocentric = - shell_p.(ab).ShellPair.monocentric && - shell_q.(cd).ShellPair.monocentric - in - *) + let () = + if debug then ( + if monocentric then + Printf.printf "Mono-centric\n" + else + Printf.printf "Multi-centric\n" + ) + in (* Compute the integral class from the primitive shell quartet *) class_indices |> Array.iteri (fun i key -> @@ -418,16 +420,19 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q | _ -> assert false in try - (* if monocentric then begin - if ( ((1 land (a.(0) + a.(3) + a.(6) + a.( 9)))=1) || - ((1 land (a.(1) + a.(4) + a.(7) + a.(10)))=1) || - ((1 land (a.(2) + a.(5) + a.(8) + a.(11)))=1) + let ax,ay,az = angMomA + and bx,by,bz = angMomB + and cx,cy,cz = angMomC + and dx,dy,dz = angMomD + in + if ( ((1 land ax+bx+cx+dx)=1) || + ((1 land ay+by+cy+dy)=1) || + ((1 land az+bz+cz+dz)=1) ) then raise NullQuartet end; - *) (* Schwartz screening *) (* let schwartz_p = diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 07aad6f..426fd8a 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -522,24 +522,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q 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 @@ -675,11 +657,12 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let norm = - let norm_coef_scale_q = shell_q.ContractedShellPair.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 + let norm_coef_scale_q = + shell_q.ContractedShellPair.norm_coef_scale + in + Array.to_list norm_coef_scale_p + |> List.map (fun v1 -> + Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) |> Array.concat in