diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index bf6dcd8..0a7595a 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -14,9 +14,6 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d) coef_prod map = - let k = 0 in - let getk a = a.(k) in - let totAngMom_a = Angular_momentum.to_int totAngMom_a and totAngMom_b = Angular_momentum.to_int totAngMom_b and totAngMom_c = Angular_momentum.to_int totAngMom_c @@ -105,29 +102,40 @@ let hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d) am.(xyz) <- am.(xyz) - 1; cm.(xyz) <- cm.(xyz) - 1; cmm.(xyz) <- cmm.(xyz) - 2; - Array.mapi (fun k _ -> -if cm.(xyz) < 0 then 0. else - let v1 = - vrr m angMom_a cm totAngMom_a (totAngMom_c-1) - and v2 = - vrr (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) - in - -. v1.(k) *. expo_d.(k) *. expo_inv_q.(k) *. (Coordinate.coord center_cd.(k) xyz) - -. v2.(k) *. (expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz)) - +. (if cmm.(xyz) < 0 then 0. else - let v1 = - vrr m angMom_a cmm totAngMom_a (totAngMom_c-2) - and v2 = - vrr (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) - in - ((float_of_int cm.(xyz)) *. expo_inv_q.(k) *. 0.5 ) *.(v1.(k) +. expo_inv_q.(k) *. v2.(k) ) ) - -. (if am.(xyz) lor cm.(xyz) < 0 then 0. else - let v = - vrr (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) - in - (float_of_int angMom_a.(xyz)) *. expo_inv_p *. expo_inv_q.(k) *. 0.5 *.v.(k) ) - ) coef_prod - ) + if cm.(xyz) < 0 then + Array.map (fun _ -> 0.) coef_prod + else + let v1 = + vrr m angMom_a cm totAngMom_a (totAngMom_c-1) + and v2 = + vrr (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) + in + let p1 = + Array.mapi (fun k _ -> + -. v1.(k) *. expo_d.(k) *. expo_inv_q.(k) *. (Coordinate.coord center_cd.(k) xyz) + -. v2.(k) *. (expo_inv_q.(k) *. (Coordinate.coord center_pq.(k) xyz)) + ) coef_prod + in + let p2 = + if cmm.(xyz) < 0 then p1 else + let v1 = + vrr m angMom_a cmm totAngMom_a (totAngMom_c-2) + and v2 = + vrr (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) + in + Array.mapi (fun k _ -> p1.(k) +. + ((float_of_int cm.(xyz)) *. expo_inv_q.(k) *. 0.5 ) + *. (v1.(k) +. expo_inv_q.(k) *. v2.(k)) + ) coef_prod + in + if (am.(xyz) < 0) || (cm.(xyz) < 0) then p2 else + let v = + vrr (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) + in + Array.mapi (fun k _ -> + p2.(k) -. (float_of_int angMom_a.(xyz)) *. expo_inv_p *. expo_inv_q.(k) *. 0.5 *.v.(k) + ) coef_prod + ) in if not found then Zmap.add map.(m) key result; @@ -142,17 +150,14 @@ if cm.(xyz) < 0 then 0. else match totAngMom_b with | 0 -> vrr m angMom_a angMom_c totAngMom_a totAngMom_c - |> getk | 1 -> let xyz = if angMom_b.(0) = 1 then 0 else if angMom_b.(1) = 1 then 1 else 2 in let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in ap.(xyz) <- ap.(xyz) + 1; let v1 = vrr m ap angMom_c (totAngMom_a+1) totAngMom_c - |> getk and v2 = vrr m angMom_a angMom_c totAngMom_a totAngMom_c - |> getk - in v1 +. v2 *. (Coordinate.coord center_ab xyz) + in Array.map2 (fun v1 v2 -> v1 +. v2 *. (Coordinate.coord center_ab xyz) ) v1 v2 | _ -> let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |] @@ -164,19 +169,18 @@ if cm.(xyz) < 0 then 0. else in ap.(xyz) <- ap.(xyz) + 1; bm.(xyz) <- bm.(xyz) - 1; - if (bm.(xyz) < 0) then 0. else + if (bm.(xyz) < 0) then Array.map (fun _ -> 0.) coef_prod else let h1 = hrr0 m ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c and h2 = hrr0 m angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c - in h1 +. h2 *. (Coordinate.coord center_ab xyz) + in Array.map2 (fun h1 h2 -> h1 +. h2 *. (Coordinate.coord center_ab xyz)) h1 h2 and hrr m angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d = match (totAngMom_b, totAngMom_d) with | (0,0) -> vrr m angMom_a angMom_c totAngMom_a totAngMom_c - |> getk | (_,0) -> hrr0 m angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c | (_,_) -> let cp = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] @@ -193,8 +197,7 @@ if cm.(xyz) < 0 then 0. else hrr m angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) , hrr m angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) in - Array.map (fun center_cd -> h1 +. h2 *. (Coordinate.coord center_cd xyz)) center_cd - |> getk + Array.mapi (fun k center_cd -> h1.(k) +. h2.(k) *. (Coordinate.coord center_cd xyz)) center_cd in hrr m angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d @@ -302,17 +305,18 @@ Array.iter (fun shell_ab -> center_pq,coef_prod) ) shell_q in - Array.iteri (fun cd shell_cd -> - try - let (zero_m_array, expo_inv, d, center_cd, - center_pq,coef_prod) = common.(cd) - in - let zero_m_array = [| zero_m_array |] - and expo_inv = [| expo_inv |] - and d = [| d |] - and center_cd = [| center_cd |] - and center_pq = [| center_pq |] - and coef_prod = [| coef_prod |] + let zero_m_array = Array.map (fun (zero_m_array, expo_inv, d, center_cd, + 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) -> expo_inv ) common + and d = Array.map (fun (zero_m_array, expo_inv, d, center_cd, + center_pq,coef_prod) -> d) common + and center_cd = Array.map (fun (zero_m_array, expo_inv, d, center_cd, + 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) -> center_pq) common + and coef_prod = Array.map (fun (zero_m_array, expo_inv, d, center_cd, + center_pq,coef_prod) -> coef_prod) common in let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in (* Compute the integral class from the primitive shell quartet *) @@ -324,26 +328,25 @@ Array.iter (fun shell_ab -> [| a.(6) ; a.(7) ; a.(8) |], [| a.(9) ; a.(10) ; a.(11) |] ) in - try - - let norm = + let norm = + Array.map (fun shell_cd -> shell_ab.Shell_pair.norm_fun angMomA angMomB *. shell_cd.Shell_pair.norm_fun angMomC angMomD - in - let integral = chop norm (fun () -> - hvrr_two_e 0 (angMomA, angMomB, angMomC, angMomD) - (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, - Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) - (maxm, zero_m_array) - (Contracted_shell.expo shell_b b, d) - (shell_ab.Shell_pair.expo_inv, expo_inv) - (shell_ab.Shell_pair.center_ab, center_cd, center_pq) - coef_prod map ) - in - contracted_class.(i) <- contracted_class.(i) +. integral - with NullQuartet -> () + ) shell_q + in + let integral = + hvrr_two_e 0 (angMomA, angMomB, angMomC, angMomD) + (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, + Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) + (maxm, zero_m_array) + (Contracted_shell.expo shell_b b, d) + (shell_ab.Shell_pair.expo_inv, expo_inv) + (shell_ab.Shell_pair.center_ab, center_cd, center_pq) + coef_prod map + |> Array.map2 (fun x y -> x *. y ) norm + in + let x = Array.fold_left (+.) 0. integral in + contracted_class.(i) <- contracted_class.(i) +. x ) class_indices - with NullQuartet -> () - ) shell_q ) shell_p end; @@ -364,4 +367,3 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = in contracted_class_shell_pairs ~zero_m shell_p shell_q -