diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index 632def2..68590c3 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -16,6 +16,7 @@ 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 *) } @@ -94,7 +95,11 @@ let create ?cutoff p_a p_b = let monocentric = Contracted_shell.center p_a = Contracted_shell.center p_b in - Some ShellPair.{ i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; monocentric } + let totAngMomInt = + (Angular_momentum.to_int (Contracted_shell.totAngMom p_a)) + + (Angular_momentum.to_int (Contracted_shell.totAngMom p_b)) + in + Some ShellPair.{ i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; monocentric ; totAngMomInt} with | Null_contribution -> None ) @@ -112,7 +117,8 @@ let create ?cutoff p_a p_b = { shell_a = p_a ; shell_b = p_b ; coef ; expo_inv ; shell_pairs ; center_ab=shell_pairs.(0).center_ab; - norm_coef_scale ; norm_sq=shell_pairs.(0).norm_sq + norm_coef_scale ; norm_sq=shell_pairs.(0).norm_sq; + totAngMomInt = shell_pairs.(0).ShellPair.totAngMomInt; } diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 58ee371..0f9e24b 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -152,18 +152,18 @@ let of_basis basis = in let swap = - Array.length sp > Array.length sq + Array.length sp > Array.length sq in (* Compute all the integrals of the class *) let cls = if swap then - if Array.length sp < 2 then + if (Array.length sp) + (Array.length sq) = 2 then contracted_class_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p else contracted_class_shell_pairs_vec ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p else - if Array.length sq < 2 then + if (Array.length sp) + (Array.length sq) = 2 then contracted_class_shell_pairs ~schwartz_p ~schwartz_q shell_p shell_q else contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q diff --git a/Basis/ShellPair.ml b/Basis/ShellPair.ml index cb9736e..98d5a99 100644 --- a/Basis/ShellPair.ml +++ b/Basis/ShellPair.ml @@ -11,6 +11,7 @@ type t = { norm_coef: float; (* norm_coef_a * norm_coef_b *) coef : float; (* norm_coef * coef_a * coef_b * g, with g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *) + totAngMomInt : int ; i : int; j : int; shell_a : Contracted_shell.t; diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 64fbdc7..ddb0891 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -318,9 +318,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q and sq = shell_q.ContractedShellPair.shell_pairs 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) + shell_p.ContractedShellPair.totAngMomInt + + shell_q.ContractedShellPair.totAngMomInt in (* Pre-computation of integral class indices *) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 01d45c7..375e5dc 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -21,10 +21,9 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (expo_b, expo_d) (expo_inv_p, expo_inv_q) (center_ab, center_cd, center_pq) - coef_prod map_1d map_2d np nq + map_1d map_2d np nq = - let empty = Array.make nq 0. in @@ -38,27 +37,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (** Vertical recurrence relations *) let rec vrr0_v m angMom_a = function - (* - | 1 -> - let xyz = - match angMom_a with - | (1,_,_) -> 0 - | (_,1,_) -> 1 - | _ -> 2 - in - Some ( - Array.init np (fun l -> - Array.init nq (fun k -> - let f = expo_b.(l) *. (Coordinate.coord center_ab xyz) in - coef_prod.(l).(k) *. expo_inv_p.(l) *. - (center_pq.(xyz).(l).(k) *. zero_m_array.(m+1).(l).(k) - -. f *. zero_m_array.(m).(l).(k) ) )) - ) - *) - | 0 -> Some ( - Array.init np (fun l -> - Array.init nq (fun k -> - zero_m_array.(m).(l).(k) *. coef_prod.(l).(k) ) ) ) + | 0 -> Some zero_m_array.(m) | totAngMom_a -> let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in @@ -97,14 +76,12 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in match v1_top with | Some v1_top -> - v1_top.(l) - |> Array.map (fun x -> f *. x) + v1_top.(l) |> Array.map (fun x -> f *. x) | None -> empty in let p1 = match p1_top with - | Some p1_top -> - p1_top.(l) + | Some p1_top -> p1_top.(l) | _ -> assert false in let p1 = @@ -131,12 +108,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in Zmap.add map_1d.(m) key result; result - and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c = match (totAngMom_a, totAngMom_c) with - | (i,0) -> - vrr0_v m angMom_a totAngMom_a + | (i,0) -> vrr0_v m angMom_a totAngMom_a | (_,_) -> let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c)) @@ -145,7 +120,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) try Zmap.find map_2d.(m) key with | Not_found -> let result = - begin + begin let am, cm, cmm, axyz, cxyz, xyz = let (aax, aay, aaz) = angMom_a and (acx, acy, acz) = angMom_c in @@ -284,26 +259,14 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | 0 -> begin match (totAngMom_a, totAngMom_c) with - | (0,0) -> - begin - let result = ref 0. in - for l=0 to np-1 do - for k=0 to nq-1 do - result := !result +. zero_m_array.(0).(l).(k) *. coef_prod.(l).(k) - done - done; - !result - end - | (_,0) -> - begin - match vrr0_v 0 angMom_a totAngMom_a with - | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix - | None -> 0. - end + | (0,0) -> Array.fold_left (fun accu c -> + accu +. Array.fold_left (+.) 0. c) 0. zero_m_array.(0) | (_,_) -> - match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with - | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix - | None -> 0. + begin + match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with + | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix + | None -> 0. + end end | 1 -> let (aax, aay, aaz) = angMom_a in @@ -376,13 +339,16 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in let h1 = hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) - and h2 = - hrr_v angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) in let f = (Coordinate.coord center_cd xyz) in - h1 +. f *. h2 + if (abs_float f < cutoff) then + h1 + else + let h2 = + hrr_v angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) + in h1 +. f *. h2 in - hrr_v + hrr_v (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)) @@ -405,9 +371,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q and sq = shell_q.ContractedShellPair.shell_pairs 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) + shell_p.ContractedShellPair.totAngMomInt + + shell_q.ContractedShellPair.totAngMomInt in (* Pre-computation of integral class indices *) @@ -543,7 +508,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (* Transpose result *) for m=0 to maxm do for cd=0 to nq-1 do - result.(m).(ab).(cd) <- zero_m_array_tmp.(cd).(m) + result.(m).(ab).(cd) <- zero_m_array_tmp.(cd).(m) *. coef.(ab).(cd) done done ) sp; @@ -559,8 +524,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q |> Array.concat in - let map_1d = Array.init maxm (fun _ -> Zmap.create (4*maxm)) - and map_2d = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) + let map_1d = Array.init (maxm+1) (fun _ -> Zmap.create (4*maxm)) + and map_2d = Array.init (maxm+1) (fun _ -> Zmap.create (Array.length class_indices)) in (* Compute the integral class from the primitive shell quartet *) Array.iteri (fun i key -> @@ -580,7 +545,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (expo_inv_p, expo_inv_q) (shell_p.ContractedShellPair.center_ab, shell_q.ContractedShellPair.center_ab, center_pq) - coef map_1d map_2d np nq + map_1d map_2d np nq in contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i) ) class_indices