diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 46d5641..ae642d0 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -15,19 +15,23 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (expo_b, expo_d) (expo_inv_p, expo_inv_q) (center_ab, center_cd, center_pq) - map_1d map_2d map_1d' map_2d' + (center_pa, center_qc) + map_1d map_2d = (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *) + (* if angMom_a.tot + angMom_b.tot < angMom_c.tot + angMom_d.tot then hvrr_two_e (angMom_c, angMom_d, angMom_a, angMom_b) zero_m_array (expo_d, expo_b) (expo_inv_q, expo_inv_p) (center_cd, center_ab, (Coordinate.neg center_pq) ) - map_1d' map_2d' map_1d map_2d + (center_qc, center_pa) + map_1d map_2d else + *) let maxm = angMom_a.tot + angMom_b.tot + angMom_c.tot + angMom_d.tot in let maxsze = maxm+1 in @@ -98,9 +102,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) try Zmap.find map_2d key with | Not_found -> let result = - (* Invariant : angMom_c.tot > 0 so cm.tot >= 0 *) - let xyz = get_xyz angMom_c in - let cm = Powers.decr xyz angMom_c in + (* angMom_c.tot > 0 so cm.tot >= 0 *) + let xyz = get_xyz angMom_c in + let cm = Powers.decr xyz angMom_c in let cmxyz = Powers.get xyz cm in let axyz = Powers.get xyz angMom_a in @@ -162,19 +166,65 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and trr angMom_a angMom_c = match (angMom_a.tot, angMom_c.tot) with - | (i,0) -> if (i>0) then vrr0 angMom_a - else zero_m_array + | (i,0) -> if (i>0) then (vrr0 angMom_a).(0) + else zero_m_array.(0) | (_,_) -> + (* let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c)) in try Zmap.find map_2d key with | Not_found -> + *) let result = - let xyz = get_xyz angMom_c in - let cm = Powers.decr xyz angMom_c in - let cmxyz = Powers.get xyz cm in - let axyz = Powers.get xyz angMom_a in - *) + let xyz = get_xyz angMom_c in + let axyz = Powers.get xyz angMom_a in + let cm = Powers.decr xyz angMom_c in + let cmxyz = Powers.get xyz cm in + + let expo_inv_q_over_p = expo_inv_q /. expo_inv_p in + let f = + Coordinate.get xyz center_qc +. expo_inv_q_over_p *. + (Coordinate.get xyz center_pa) + in + let result = + if abs_float f < cutoff then 0. else + let v1 = trr angMom_a cm in + f *. v1 + in + let result = + if axyz < 1 then result else + let f = 0.5 *. (float_of_int axyz) *. expo_inv_q in + if abs_float f < cutoff then result else + let am = Powers.decr xyz angMom_a in + let v2 = trr am cm in + result +. f *. v2 + in + let result = + if cmxyz < 1 then result else + let f = 0.5 *. (float_of_int cmxyz) *. expo_inv_q in + if abs_float f < cutoff then 0. else + let cmm = Powers.decr xyz cm in + let v3 = trr angMom_a cmm in + result +. f *. v3 + in + let result = + if cmxyz < 0 then result else + let f = -. expo_inv_q_over_p in + let ap = Powers.incr xyz angMom_a in + let v4 = trr ap cm in + result +. v4 *. f + in + result + in + (* + Zmap.add map_2d key result; + *) + result + *) + + + + @@ -183,7 +233,6 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and hrr0 angMom_a angMom_b angMom_c = match angMom_b.tot with - | 0 -> (vrr angMom_a angMom_c).(0) | 1 -> let xyz = get_xyz angMom_b in let ap = Powers.incr xyz angMom_a in @@ -192,6 +241,10 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) if (abs_float f2 < cutoff) then v1.(0) else let v2 = vrr angMom_a angMom_c in v1.(0) +. f2 *. v2.(0) + | 0 -> (vrr angMom_a angMom_c).(0) + (* + | 0 -> trr angMom_a angMom_c + *) | _ -> let xyz = get_xyz angMom_b in let bxyz = Powers.get xyz angMom_b in @@ -212,6 +265,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) | (_,0) -> if (angMom_b.tot = 0) then (vrr angMom_a angMom_c).(0) + (* + trr angMom_a angMom_c + *) else hrr0 angMom_a angMom_b angMom_c | (_,_) -> @@ -306,8 +362,6 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let d = shell_q.ContractedShellPair.shell_pairs.(cd).ShellPair.j in let map_1d = Zmap.create (4*maxm) in let map_2d = Zmap.create (Array.length class_indices) in - let map_1d' = Zmap.create (4*maxm) in - let map_2d' = Zmap.create (Array.length class_indices) in let norm_coef_scale = Array.to_list norm_coef_scale_p |> List.map (fun v1 -> @@ -361,6 +415,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let norm = norm_coef_scale.(i) in let coef_prod = coef_prod *. norm in + let integral = hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) zero_m_array @@ -368,7 +423,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (shell_p.ContractedShellPair.expo_inv.(ab), shell_q.ContractedShellPair.expo_inv.(cd) ) (sp.(ab).ShellPair.center_ab, sq.(cd).ShellPair.center_ab, center_pq) - map_1d map_2d map_1d' map_2d' + (sp.(ab).ShellPair.center_a , sq.(cd).ShellPair.center_a) + map_1d map_2d in contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral with NullQuartet -> ()