diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 36ecd9b..d275ed9 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -70,9 +70,9 @@ let index i j k l = let of_basis basis = - let to_int_tuple x = + let to_powers x = let open Zkey in - match to_int_tuple Kind_3 x with + match to_powers Kind_3 x with | Three x -> x | _ -> assert false in @@ -173,21 +173,21 @@ let of_basis basis = (* Write the data in the output file *) Array.iteri (fun i_c powers_i -> let i_c = (Contracted_shell.index shell.(i)) + i_c + 1 in - let xi = to_int_tuple powers_i in + let xi = to_powers powers_i in Array.iteri (fun j_c powers_j -> let j_c = (Contracted_shell.index shell.(j)) + j_c + 1 in - let xj = to_int_tuple powers_j in + let xj = to_powers powers_j in Array.iteri (fun k_c powers_k -> let k_c = (Contracted_shell.index shell.(k)) + k_c + 1 in - let xk = to_int_tuple powers_k in + let xk = to_powers powers_k in Array.iteri (fun l_c powers_l -> let l_c = (Contracted_shell.index shell.(l)) + l_c + 1 in - let xl = to_int_tuple powers_l in + let xl = to_powers powers_l in let key = if swap then - Zkey.of_int_tuple (Zkey.Twelve (xk,xl,xi,xj)) + Zkey.of_powers (Zkey.Twelve (xk,xl,xi,xj)) else - Zkey.of_int_tuple (Zkey.Twelve (xi,xj,xk,xl)) + Zkey.of_powers (Zkey.Twelve (xi,xj,xk,xl)) in let value = Zmap.find cls key diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index 526ea36..e10a1ae 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -1,6 +1,7 @@ open Util open Constants open Lacaml.D +open Powers open Coordinate type t = Mat.t @@ -108,9 +109,9 @@ let contracted_class shell_a shell_b : float Zmap.t = (** Create kinetic energy matrix *) let of_basis basis = - let to_int_tuple x = + let to_powers x = let open Zkey in - match to_int_tuple Kind_3 x with + match to_powers Kind_3 x with | Three x -> x | _ -> assert false in @@ -131,12 +132,12 @@ let of_basis basis = Array.iteri (fun j_c powers_j -> let j_c = Contracted_shell.index shell.(j) + j_c + 1 in - let xj = to_int_tuple powers_j in + let xj = to_powers powers_j in Array.iteri (fun i_c powers_i -> let i_c = Contracted_shell.index shell.(i) + i_c + 1 in - let xi = to_int_tuple powers_i in + let xi = to_powers powers_i in let key = - Zkey.of_int_tuple (Zkey.Six (xi,xj)) + Zkey.of_powers (Zkey.Six (xi,xj)) in let value = try Zmap.find cls key diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml index d28d9bf..cd85b00 100644 --- a/Basis/NucInt.ml +++ b/Basis/NucInt.ml @@ -36,9 +36,9 @@ exception NullIntegral let of_basis_nuclei basis nuclei = - let to_int_tuple x = + let to_powers x = let open Zkey in - match to_int_tuple Kind_3 x with + match to_powers Kind_3 x with | Three x -> x | _ -> assert false in @@ -72,12 +72,12 @@ let of_basis_nuclei basis nuclei = (* Write the data in the output file *) Array.iteri (fun i_c powers_i -> let i_c = (Contracted_shell.index shell.(i)) + i_c + 1 in - let xi = to_int_tuple powers_i in + let xi = to_powers powers_i in Array.iteri (fun j_c powers_j -> let j_c = (Contracted_shell.index shell.(j)) + j_c + 1 in - let xj = to_int_tuple powers_j in + let xj = to_powers powers_j in let key = - Zkey.of_int_tuple (Zkey.Six (xi,xj)) + Zkey.of_powers (Zkey.Six (xi,xj)) in let value = Zmap.find cls key diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index 38f72e5..bd4c114 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -1,5 +1,6 @@ open Util open Constants +open Powers open Coordinate exception NullPair @@ -12,56 +13,57 @@ let chop f g = (** Horizontal and Vertical Recurrence Relations (HVRR) *) -let hvrr_one_e - (angMom_a, angMom_b) (totAngMom_a, totAngMom_b) - (maxm, zero_m_array) (expo_b) (expo_inv_p) (center_ab, center_pa, center_pc) +let hvrr_one_e (angMom_a, angMom_b) + zero_m_array (expo_b) (expo_inv_p) (center_ab, center_pa, center_pc) map = - let totAngMom_a = Angular_momentum.to_int totAngMom_a - and totAngMom_b = Angular_momentum.to_int totAngMom_b - in - let maxm = totAngMom_a+totAngMom_b in + let maxm = angMom_a.tot + angMom_b.tot in let maxsze = maxm+1 in let empty = Array.make maxsze 0. in + let get_xyz angMom = + match angMom with + | { y=0 ; z=0 ; _ } -> X + | { z=0 ; _ } -> Y + | _ -> Z + in + + (** Vertical recurrence relations *) - let rec vrr angMom_a totAngMom_a = - let ax,ay,az = angMom_a in - if ax < 0 || ay < 0 || az < 0 then - empty + let rec vrr angMom_a = + let { x=ax ; y=ay ; z=az } = angMom_a in + if ax < 0 || ay < 0 || az < 0 then raise Exit else - match totAngMom_a with + match angMom_a.tot with | 0 -> zero_m_array | _ -> - let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in + let key = Zkey.of_powers (Zkey.Three angMom_a) in try Zmap.find map key with | Not_found -> let result = - let am, amm, amxyz, xyz = - match angMom_a with - | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, X - | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, Y - | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, Z - in + let xyz = get_xyz angMom_a in + let am = Powers.decr xyz angMom_a in + let amxyz = Powers.get xyz am in if amxyz < 0 then empty else let f1 = Coordinate.get xyz center_pa and f2 = expo_inv_p *. (Coordinate.get xyz center_pc) in if amxyz < 1 then let v1 = - vrr am (totAngMom_a-1) + vrr am in Array.init maxsze (fun m -> if m = maxm then (f1 *. v1.(m) ) else (f1 *. v1.(m) ) -. f2 *. v1.(m+1) ) else let v3 = - vrr amm (totAngMom_a-2) + let amm = Powers.decr xyz am in + vrr amm in let v1 = - vrr am (totAngMom_a-1) + vrr am in let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in Array.init maxsze (fun m -> f1 *. v1.(m) -. @@ -75,43 +77,33 @@ let hvrr_one_e (** Horizontal recurrence relations *) - and hrr angMom_a angMom_b totAngMom_a totAngMom_b = + and hrr angMom_a angMom_b = - let bx,by,bz = angMom_b in - if bx < 0 || by < 0 || bz < 0 then 0. + let { x=bx ; y=by ; z=bz } = angMom_b in + if bx < 0 || by < 0 || bz < 0 then raise Exit else - match totAngMom_b with - | 0 -> (vrr angMom_a totAngMom_a).(0) + match angMom_b.tot with + | 0 -> (vrr angMom_a).(0) | _ -> - let angMom_ax, angMom_ay, angMom_az = angMom_a - and angMom_bx, angMom_by, angMom_bz = angMom_b in - let bxyz, xyz = - match angMom_b with - | (_,0,0) -> angMom_bx, X - | (_,_,0) -> angMom_by, Y - | (_,_,_) -> angMom_bz, Z - in + let xyz = get_xyz angMom_b in + let bxyz = Powers.get xyz angMom_b in if (bxyz < 1) then 0. else - let ap, bm = - match xyz with - | X -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz) - | Y -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz) - | Z -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1) - in + let ap = Powers.incr xyz angMom_a in + let bm = Powers.decr xyz angMom_b in let h1 = - hrr ap bm (totAngMom_a+1) (totAngMom_b-1) + hrr ap bm in let f2 = Coordinate.get xyz center_ab in if abs_float f2 < cutoff then h1 else let h2 = - hrr angMom_a bm totAngMom_a (totAngMom_b-1) + hrr angMom_a bm in h1 +. f2 *. h2 in - hrr angMom_a angMom_b totAngMom_a totAngMom_b + hrr angMom_a angMom_b @@ -199,7 +191,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = class_indices |> Array.iteri (fun i key -> let (angMomA,angMomB) = - match Zkey.to_int_tuple ~kind:Zkey.Kind_6 key with + match Zkey.to_powers ~kind:Zkey.Kind_6 key with | Zkey.Six x -> x | _ -> assert false in @@ -207,8 +199,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let coef_prod = coef_prod *. norm in let integral = hvrr_one_e (angMomA, angMomB) - (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b) - (maxm, zero_m_array) + zero_m_array (Contracted_shell.expo shell_b b) (shell_p.ContractedShellPair.expo_inv.(ab)) (center_ab, center_pa, center_pc) diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index e40b4c3..22ef83b 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -82,9 +82,9 @@ let contracted_class shell_a shell_b : float Zmap.t = (** Create overlap matrix *) let of_basis basis = - let to_int_tuple x = + let to_powers x = let open Zkey in - match to_int_tuple Kind_3 x with + match to_powers Kind_3 x with | Three x -> x | _ -> assert false in @@ -105,12 +105,12 @@ let of_basis basis = Array.iteri (fun j_c powers_j -> let j_c = Contracted_shell.index shell.(j) + j_c + 1 in - let xj = to_int_tuple powers_j in + let xj = to_powers powers_j in Array.iteri (fun i_c powers_i -> let i_c = Contracted_shell.index shell.(i) + i_c + 1 in - let xi = to_int_tuple powers_i in + let xi = to_powers powers_i in let key = - Zkey.of_int_tuple (Zkey.Six (xi,xj)) + Zkey.of_powers (Zkey.Six (xi,xj)) in let value = try Zmap.find cls key diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 9c85841..2abedb3 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -1,5 +1,6 @@ open Util open Constants +open Powers open Coordinate let debug=false @@ -10,34 +11,37 @@ exception NullQuartet (** Horizontal and Vertical Recurrence Relations (HVRR) *) let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) - (totAngMom_a_in, totAngMom_b_in, totAngMom_c_in, totAngMom_d_in) - (maxm, zero_m_array) + zero_m_array (expo_b, expo_d) (expo_inv_p, expo_inv_q) (center_ab, center_cd, center_pq) map_1d map_2d map_1d' map_2d' = - let maxsze = maxm+1 in - let totAngMom_a = Angular_momentum.to_int totAngMom_a_in - and totAngMom_b = Angular_momentum.to_int totAngMom_b_in - and totAngMom_c = Angular_momentum.to_int totAngMom_c_in - and totAngMom_d = Angular_momentum.to_int totAngMom_d_in - in - (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *) - if (totAngMom_a+totAngMom_b < totAngMom_c+totAngMom_d) then + 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) - (totAngMom_c_in, totAngMom_d_in, totAngMom_a_in, totAngMom_b_in) - (maxm, zero_m_array) + 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 + else - let maxm = totAngMom_a + totAngMom_b + totAngMom_c + totAngMom_d in - let empty = Array.make (maxm+1) 0. + + let maxm = angMom_a.tot + angMom_b.tot + angMom_c.tot + angMom_d.tot in + let maxsze = maxm+1 in + + let empty = Array.make (maxm+1) 0. in + + let get_xyz angMom = + match angMom with + | { y=0 ; z=0 ; _ } -> X + | { z=0 ; _ } -> Y + | _ -> Z in + + (* if debug then begin Printf.printf "\n---- %d %d %d %d ----\n" totAngMom_a totAngMom_b totAngMom_c totAngMom_d; let (x,y,z) = angMom_a in Printf.printf "%d %d %d\n" x y z; @@ -50,29 +54,30 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (get X center_cd) (get Y center_cd) (get Z center_cd) (get X center_pq) (get Y center_pq) (get Z center_pq) end; + *) (** Vertical recurrence relations *) - let rec vrr0 angMom_a totAngMom_a = + let rec vrr0 angMom_a = + (* if debug then begin let (x,y,z) = angMom_a in - Printf.printf "vrr0: %d : %d %d %d\n" totAngMom_a x y z + Printf.printf "vrr0: %d : %d %d %d\n" angMom_a.tot x y z end; + *) - match totAngMom_a with + match angMom_a.tot with | 0 -> zero_m_array | _ -> - let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in + let key = Zkey.of_powers (Zkey.Three angMom_a) in try Zmap.find map_1d key with | Not_found -> let result = - let am, amm, amxyz, xyz = - match angMom_a with - | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, X - | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, Y - | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, Z - in + let xyz = get_xyz angMom_a in + let am = Powers.decr xyz angMom_a in + let amxyz = Powers.get xyz am in + if amxyz < 0 then empty else let f1 = expo_inv_p *. (Coordinate.get xyz center_pq) and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab) @@ -81,7 +86,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) if amxyz < 1 then begin let v1 = - vrr0 am (totAngMom_a-1) + vrr0 am in for m=0 to maxm-1 do result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) @@ -90,12 +95,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) end else begin - let v3 = - vrr0 amm (totAngMom_a-2) - in - let v1 = - vrr0 am (totAngMom_a-1) - in + let amm = Powers.decr xyz am in + let v3 = vrr0 amm in + let v1 = vrr0 am in let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in for m=0 to maxm-1 do result.(m) <- f1 *. v1.(m+1) -. f2 *. v1.(m) @@ -108,51 +110,38 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) result - and vrr angMom_a angMom_c totAngMom_a totAngMom_c = + and vrr angMom_a angMom_c = + (* if debug then begin let angMom_ax, angMom_ay, angMom_az = angMom_a in let angMom_cx, angMom_cy, angMom_cz = angMom_c in - Printf.printf "vrr : %d %d : %d %d %d %d %d %d\n" totAngMom_a totAngMom_c + Printf.printf "vrr : %d %d : %d %d %d %d %d %d\n" angMom_a.tot angMom_c.tot angMom_ax angMom_ay angMom_az angMom_cx angMom_cy angMom_cz end; + *) - match (totAngMom_a, totAngMom_c) with + match (angMom_a.tot, angMom_c.tot) with | (i,0) -> if (i>0) then - vrr0 angMom_a totAngMom_a + vrr0 angMom_a (* - OneElectronRR.hvrr_one_e (angMom_a, angMom_b) (totAngMom_a_in, totAngMom_b_in) + OneElectronRR.hvrr_one_e (angMom_a, angMom_b) (angMom_a.tot, angMom_b.tot) (maxm, zero_m_array) (expo_b) (expo_inv_p) (center_ab, center_pq, center_ab) map_1d *) else zero_m_array | (_,_) -> - let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c) ) in + let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c)) in try Zmap.find map_2d key with | Not_found -> let result = - let am, cm, cmm, axyz, cmxyz, xyz = - 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) -> - (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, X - | (_,_,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, Y - | _ -> - (angMom_ax, angMom_ay, angMom_az-1), - (angMom_cx, angMom_cy, angMom_cz-1), - (angMom_cx, angMom_cy, angMom_cz-2), - angMom_az,angMom_cz-1, Z - in + 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 + if cmxyz < 0 then empty else let f1 = @@ -164,7 +153,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then begin let v1 = - vrr angMom_a cm totAngMom_a (totAngMom_c-1) + vrr angMom_a cm in for m=0 to maxm-1 do result.(m) <- f1 *. v1.(m) -. f2 *. v1.(m+1) ; @@ -180,7 +169,8 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (abs_float (f3 *. expo_inv_q) > cutoff) then begin let v3 = - vrr angMom_a cmm totAngMom_a (totAngMom_c-2) + let cmm = Powers.decr xyz cm in + vrr angMom_a cmm in for m=0 to maxm-1 do result.(m) <- result.(m) +. @@ -191,12 +181,13 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) end; if (axyz > 0) && (cmxyz >= 0) then begin + let am = Powers.decr xyz angMom_a in let f5 = (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 in if (abs_float f5 > cutoff) then let v5 = - vrr am cm (totAngMom_a-1) (totAngMom_c-1) + vrr am cm in for m=0 to maxm-1 do result.(m) <- result.(m) -. f5 *. v5.(m+1) @@ -211,8 +202,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (** Horizontal recurrence relations *) - and hrr0 angMom_a angMom_b angMom_c - totAngMom_a totAngMom_b totAngMom_c = + and hrr0 angMom_a angMom_b angMom_c = (* if debug then @@ -221,93 +211,70 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and angMom_bx, angMom_by, angMom_bz = angMom_b and angMom_cx, angMom_cy, angMom_cz = angMom_c in Printf.printf "hrr0: %d %d %d : %d %d %d %d %d %d %d %d %d\n" - totAngMom_a totAngMom_b totAngMom_c angMom_ax angMom_ay angMom_az angMom_bx angMom_by angMom_bz angMom_cx angMom_cy angMom_cz end; *) - match totAngMom_b with - | 0 -> (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0) + match angMom_b.tot with + | 0 -> (vrr angMom_a angMom_c).(0) | 1 -> - let angMom_ax, angMom_ay, angMom_az = angMom_a in - let ap, xyz = - match angMom_b with - | (1,_,_) -> (angMom_ax+1,angMom_ay,angMom_az), X - | (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), Y - | _ -> (angMom_ax,angMom_ay,angMom_az+1), Z - in + let xyz = get_xyz angMom_b in + let ap = Powers.incr xyz angMom_a in let v1 = - vrr ap angMom_c (totAngMom_a+1) totAngMom_c + vrr ap angMom_c in let f2 = (Coordinate.get xyz center_ab) in if (abs_float f2 < cutoff) then v1.(0) else let v2 = - vrr angMom_a angMom_c totAngMom_a totAngMom_c + vrr angMom_a angMom_c in v1.(0) +. f2 *. v2.(0) | _ -> - let angMom_ax, angMom_ay, angMom_az = angMom_a - and angMom_bx, angMom_by, angMom_bz = angMom_b in - let bxyz, xyz = - match angMom_b with - | (_,0,0) -> angMom_bx, X - | (_,_,0) -> angMom_by, Y - | (_,_,_) -> angMom_bz, Z - in + let xyz = get_xyz angMom_b in + let bxyz = Powers.get xyz angMom_b in if (bxyz < 1) then 0. else - let ap, bm = - match xyz with - | X -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz) - | Y -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz) - | Z -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1) - in + let ap = Powers.incr xyz angMom_a in + let bm = Powers.decr xyz angMom_b in let h1 = - hrr0 ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c + hrr0 ap bm angMom_c in let f2 = (Coordinate.get xyz center_ab) in if (abs_float f2 < cutoff) then h1 else let h2 = - hrr0 angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c + hrr0 angMom_a bm angMom_c in h1 +. f2 *. h2 - and hrr angMom_a angMom_b angMom_c angMom_d - totAngMom_a totAngMom_b totAngMom_c totAngMom_d = + and hrr angMom_a angMom_b angMom_c angMom_d = - match (totAngMom_b, totAngMom_d) with - | (_,0) -> if (totAngMom_b = 0) then - (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0) + match (angMom_b.tot, angMom_d.tot) with + | (_,0) -> if (angMom_b.tot = 0) then + (vrr angMom_a angMom_c).(0) else - hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c + hrr0 angMom_a angMom_b angMom_c | (_,_) -> - let (angMom_cx, angMom_cy, angMom_cz) = angMom_c - and (angMom_dx, angMom_dy, angMom_dz) = angMom_d in - let cp, dm, xyz = - match angMom_d with - | (_,0,0) -> (angMom_cx+1, angMom_cy, angMom_cz), (angMom_dx-1, angMom_dy, angMom_dz), X - | (_,_,0) -> (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), Y - | _ -> (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), Z - in + let xyz = get_xyz angMom_d in + let cp = Powers.incr xyz angMom_c in + let dm = Powers.decr xyz angMom_d in let h1 = - hrr angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) + hrr angMom_a angMom_b cp dm in let f2 = Coordinate.get xyz center_cd in if (abs_float f2 < cutoff) then h1 else let h2 = - hrr angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) + hrr angMom_a angMom_b angMom_c dm in h1 +. f2 *. h2 in hrr angMom_a angMom_b angMom_c angMom_d - totAngMom_a totAngMom_b totAngMom_c totAngMom_d @@ -400,25 +367,21 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (* Compute the integral class from the primitive shell quartet *) class_indices |> Array.iteri (fun i key -> - let (angMomA,angMomB,angMomC,angMomD) = - match Zkey.to_int_tuple ~kind:Zkey.Kind_12 key with + let (angMom_a,angMom_b,angMom_c,angMom_d) = + match Zkey.to_powers ~kind:Zkey.Kind_12 key with | Zkey.Twelve x -> x | _ -> assert false in try if monocentric then begin - 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) + if ( ((1 land angMom_a.x+angMom_b.x+angMom_c.x+angMom_d.x)=1) || + ((1 land angMom_a.y+angMom_b.y+angMom_c.y+angMom_d.y)=1) || + ((1 land angMom_a.z+angMom_b.z+angMom_c.z+angMom_d.z)=1) ) then raise NullQuartet end; + (* (* Schwartz screening *) if (maxm > 2) then @@ -448,10 +411,8 @@ 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 (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) + hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) + zero_m_array (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d) (shell_p.ContractedShellPair.expo_inv.(ab), shell_q.ContractedShellPair.expo_inv.(cd) ) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 2aca64f..066316d 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -18,25 +18,13 @@ let at_least_one_valid arr = (** Horizontal and Vertical Recurrence Relations (HVRR) *) let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) - (totAngMom_a, totAngMom_b, totAngMom_c, totAngMom_d) - (maxm, zero_m_array) + zero_m_array (expo_b, expo_d) (expo_inv_p, expo_inv_q) (center_ab, center_cd, center_pq) map_1d map_2d np nq = - let angMom_a = Powers.of_int_tuple angMom_a - and angMom_b = Powers.of_int_tuple angMom_b - and angMom_c = Powers.of_int_tuple angMom_c - and angMom_d = Powers.of_int_tuple angMom_d - in - let totAngMom_a = angMom_a.tot - and totAngMom_b = angMom_b.tot - and totAngMom_c = angMom_c.tot - and totAngMom_d = angMom_d.tot - in - let get_xyz angMom = match angMom with | { y=0 ; z=0 ; _ } -> X @@ -45,126 +33,115 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in (** Vertical recurrence relations *) - let rec vrr0_v m angMom_a = function + let rec vrr0_v m angMom_a = + match angMom_a.tot with | 0 -> Some zero_m_array.(m) - | totAngMom_a -> - let key = Zkey.of_int_tuple (Zkey.Three (Powers.to_int_tuple angMom_a)) + | _ -> + let key = Zkey.of_powers (Zkey.Three angMom_a) in try Zmap.find map_1d.(m) key with | Not_found -> let result = - try - let xyz = get_xyz angMom_a in - let am = - try Powers.decr xyz angMom_a - with Invalid_argument _ -> raise Exit + let xyz = get_xyz angMom_a in + let am = Powers.decr xyz angMom_a in + let amxyz = Powers.get xyz am in + if amxyz >= 0 then + begin + let cab = Coordinate.get xyz center_ab in + let v1_top, p1_top = + if abs_float cab < cutoff then + None, + vrr0_v (m+1) am + else + vrr0_v m am, vrr0_v (m+1) am in - (* - if amxyz < 0 then - raise Exit - else - *) - begin - let cab = Coordinate.get xyz center_ab in - let v1_top, p1_top = - if abs_float cab < cutoff then - None, - vrr0_v (m+1) am (totAngMom_a-1) - else - vrr0_v m am (totAngMom_a-1), - vrr0_v (m+1) am (totAngMom_a-1) - in - let v1_top2, p1_top2 = - try - (* - if amxyz < 1 then None, None else - *) - let amm = Powers.decr xyz am in - vrr0_v m amm (totAngMom_a-2), - vrr0_v (m+1) amm (totAngMom_a-2) - with Invalid_argument _ -> None, None - in - - let result = Array.make_matrix np nq 0. in - let p0 = - match p1_top with - | Some p1_top -> p1_top - | _ -> assert false - in - begin - match v1_top with - | None -> () - | Some v0 -> - Array.iteri (fun l result_l -> - let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab - and v0_l = v0.(l) - in - Array.iteri (fun k v0_lk -> + let v1_top2, p1_top2 = + if amxyz < 1 then None, None else + let amm = Powers.decr xyz am in + vrr0_v m amm, vrr0_v (m+1) amm + in + + let result = Array.make_matrix np nq 0. in + let p0 = + match p1_top with + | Some p1_top -> p1_top + | _ -> assert false + in + begin + match v1_top with + | None -> () + | Some v0 -> + Array.iteri (fun l result_l -> + let f0 = -. expo_b.(l) *. expo_inv_p.(l) *. cab + and v0_l = v0.(l) + in + Array.iteri (fun k v0_lk -> result_l.(k) <- v0_lk *. f0) v0_l ) result - end; - let amxyz = Powers.get xyz am in - if amxyz < 1 then - Array.iteri (fun l result_l -> - let expo_inv_p_l = expo_inv_p.(l) - and center_pq_xyz_l = (center_pq xyz).(l) - and result_l = result.(l) - and p0_l = p0.(l) - in - Array.iteri (fun k p0_lk -> - result_l.(k) <- result_l.(k) - +. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk - ) p0_l ) result - else - begin - let v1 = - match v1_top2 with - | Some v1_top2 -> v1_top2 - | None -> assert false + end; + let amxyz = Powers.get xyz am in + if amxyz < 1 then + Array.iteri (fun l result_l -> + let expo_inv_p_l = expo_inv_p.(l) + and center_pq_xyz_l = (center_pq xyz).(l) + and result_l = result.(l) + and p0_l = p0.(l) + in + Array.iteri (fun k p0_lk -> + result_l.(k) <- result_l.(k) + +. expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk + ) p0_l ) result + else + begin + let v1 = + match v1_top2 with + | Some v1_top2 -> v1_top2 + | None -> assert false + in + let v2 = + match p1_top2 with + | Some p1_top2 -> p1_top2 + | None -> assert false + in + Array.iteri (fun l result_l -> + let f = float_of_int amxyz *. expo_inv_p.(l) *. 0.5 + and expo_inv_p_l = expo_inv_p.(l) + and center_pq_xyz_l = (center_pq xyz).(l) + and v1_l = v1.(l) + and v2_l = v2.(l) + and result_l = result.(l) in - let v2 = - match p1_top2 with - | Some p1_top2 -> p1_top2 - | None -> assert false - in - Array.iteri (fun l result_l -> - let f = float_of_int amxyz *. expo_inv_p.(l) *. 0.5 - and expo_inv_p_l = expo_inv_p.(l) - and center_pq_xyz_l = (center_pq xyz).(l) - and v1_l = v1.(l) - and v2_l = v2.(l) - and result_l = result.(l) - in - Array.iteri (fun k p0_lk -> + Array.iteri (fun k p0_lk -> result_l.(k) <- result_l.(k) +. - expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +. - f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l) + expo_inv_p_l *. center_pq_xyz_l.(k) *. p0_lk +. + f *. (v1_l.(k) +. v2_l.(k) *. expo_inv_p_l) ) p0.(l) ) result - end; - Some result - end - with Exit -> None - in - Zmap.add map_1d.(m) key result; - result + end; + Some result + end + else + None + in + Zmap.add map_1d.(m) key result; + result - and vrr_v m angMom_a angMom_c totAngMom_a totAngMom_c = + and vrr_v m angMom_a angMom_c = - match (totAngMom_a, totAngMom_c) with - | (i,0) -> vrr0_v m angMom_a totAngMom_a + match (angMom_a.tot, angMom_c.tot) with + | (i,0) -> vrr0_v m angMom_a | (_,_) -> - let key = Zkey.of_int_tuple (Zkey.Six Powers.(to_int_tuple angMom_a, to_int_tuple angMom_c)) + let key = Zkey.of_powers (Zkey.Six (angMom_a, angMom_c)) in try Zmap.find map_2d.(m) key with | Not_found -> let result = - begin + begin let xyz = get_xyz angMom_c in - let cm = Powers.decr xyz angMom_c in - let axyz = Powers.get xyz angMom_a in + let cm = Powers.decr xyz angMom_c in + let axyz = Powers.get xyz angMom_a in let do_compute = ref false in let v1 = @@ -172,56 +149,56 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let f1 = Array.init nq (fun k -> - let x = expo_d.(k) *. expo_inv_q.(k) *. f in - if ( (not !do_compute) && (abs_float x > cutoff) ) then - do_compute := true; - x) + let x = expo_d.(k) *. expo_inv_q.(k) *. f in + if ( (not !do_compute) && (abs_float x > cutoff) ) then + do_compute := true; + x) in if (!do_compute) then - match vrr_v m angMom_a cm totAngMom_a (totAngMom_c-1) with + match vrr_v m angMom_a cm with | None -> None | Some v1 -> - begin - let result = Array.make_matrix np nq 0. in - for l=0 to np-1 do - for k=0 to nq-1 do - result.(l).(k) <- v1.(l).(k) *. f1.(k) - done - done; - Some ( + begin + let result = Array.make_matrix np nq 0. in + for l=0 to np-1 do + for k=0 to nq-1 do + result.(l).(k) <- v1.(l).(k) *. f1.(k) + done + done; + Some ( Array.init np (fun l -> - let v1_l = v1.(l) in - Array.init nq (fun k -> v1_l.(k) *. f1.(k)) - )) - end + let v1_l = v1.(l) in + Array.init nq (fun k -> v1_l.(k) *. f1.(k)) + )) + end else None in let v2 = let f2 = Array.init np (fun l -> - let cpq_l = (center_pq xyz).(l) in - Array.init nq (fun k -> - let x = expo_inv_q.(k) *. cpq_l.(k) in - if (!do_compute) then x - else (if abs_float x > cutoff then do_compute := true ; x) - ) ) + let cpq_l = (center_pq xyz).(l) in + Array.init nq (fun k -> + let x = expo_inv_q.(k) *. cpq_l.(k) in + if (!do_compute) then x + else (if abs_float x > cutoff then do_compute := true ; x) + ) ) in if (!do_compute) then - match vrr_v (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) with + match vrr_v (m+1) angMom_a cm with | None -> None | Some v2 -> - begin - for l=0 to np-1 do - let f2_l = f2.(l) - and v2_l = v2.(l) - in - for k=0 to nq-1 do - f2_l.(k) <- -. v2_l.(k) *. f2_l.(k) - done - done; - Some f2 - end + begin + for l=0 to np-1 do + let f2_l = f2.(l) + and v2_l = v2.(l) + in + for k=0 to nq-1 do + f2_l.(k) <- -. v2_l.(k) *. f2_l.(k) + done + done; + Some f2 + end else None in @@ -232,81 +209,75 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | None, Some v2 -> Some v2 | Some v1, None -> Some v1 | Some v1, Some v2 -> - begin - for l=0 to np-1 do - let v1_l = v1.(l) - and v2_l = v2.(l) - in - for k=0 to nq-1 do - v2_l.(k) <- v2_l.(k) +. v1_l.(k) - done - done; - Some v2 - end + begin + for l=0 to np-1 do + let v1_l = v1.(l) + and v2_l = v2.(l) + in + for k=0 to nq-1 do + v2_l.(k) <- v2_l.(k) +. v1_l.(k) + done + done; + Some v2 + end in let cxyz = Powers.get xyz angMom_c in let p2 = - try - let cmm = - try Powers.decr xyz cm - with Invalid_argument _ -> raise Exit - in - if Powers.get xyz cmm < 0 then - raise Exit; - + if cxyz < 2 then p1 else + let cmm = Powers.decr xyz cm in let fcm = (float_of_int (cxyz-1)) *. 0.5 in let f1 = Array.init nq (fun k -> - let x = fcm *. expo_inv_q.(k) in - if (!do_compute) then x - else (if abs_float x > cutoff then do_compute := true ; x) - ) + let x = fcm *. expo_inv_q.(k) in + if (!do_compute) then x + else (if abs_float x > cutoff then do_compute := true ; x) + ) in let v1 = if (!do_compute) then - match vrr_v m angMom_a cmm totAngMom_a (totAngMom_c-2) with + match vrr_v m angMom_a cmm with | None -> None | Some v1 -> - begin - let result = Array.make_matrix np nq 0. in - for l=0 to np-1 do - let v1_l = v1.(l) - and result_l = result.(l) - in - for k=0 to nq-1 do - result_l.(k) <- v1_l.(k) *. f1.(k) - done; + begin + let result = Array.make_matrix np nq 0. in + for l=0 to np-1 do + let v1_l = v1.(l) + and result_l = result.(l) + in + for k=0 to nq-1 do + result_l.(k) <- v1_l.(k) *. f1.(k) done; - Some result - end + done; + Some result + end else None in let v3 = let f2 = Array.init nq (fun k -> - let x = expo_inv_q.(k) *. f1.(k) in - if (!do_compute) then x - else (if abs_float x > cutoff then do_compute := true ; x) - ) + let x = expo_inv_q.(k) *. f1.(k) in + if (!do_compute) then x + else (if abs_float x > cutoff then do_compute := true ; x) + ) in if (!do_compute) then - match vrr_v (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) with + match vrr_v (m+1) angMom_a cmm with | None -> None | Some v3 -> - begin - let result = Array.make_matrix np nq 0. in - for l=0 to np-1 do - let v3_l = v3.(l) - and result_l = result.(l) - in - for k=0 to nq-1 do - result_l.(k) <- v3_l.(k) *. f2.(k) - done - done; - Some result - end + begin + let result = Array.make_matrix np nq 0. in + for l=0 to np-1 do + let v3_l = v3.(l) + and result_l = result.(l) + in + for k=0 to nq-1 do + result_l.(k) <- v3_l.(k) *. f2.(k) + done + done; + Some result + end else None in match p1, v1, v3 with @@ -315,88 +286,81 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | None, Some v1, None -> Some v1 | None, None, Some v3 -> Some v3 | Some p1, Some v1, Some v3 -> - begin - for l=0 to np-1 do - let v3_l = v3.(l) - and v1_l = v1.(l) - and p1_l = p1.(l) - in - for k=0 to nq-1 do - v3_l.(k) <- p1_l.(k) +. v1_l.(k) +. v3_l.(k) - done - done; - Some v3 - end - | Some p1, Some v1, None -> - begin - for l=0 to np-1 do - let v1_l = v1.(l) - and p1_l = p1.(l) - in - for k=0 to nq-1 do - p1_l.(k) <- v1_l.(k) +. p1_l.(k) - done - done; - Some p1 - end - | Some p1, None, Some v3 -> - begin - for l=0 to np-1 do - let v3_l = v3.(l) - and p1_l = p1.(l) - in - for k=0 to nq-1 do - p1_l.(k) <- p1_l.(k) +. v3_l.(k) - done - done; - Some p1 - end - | None , Some v1, Some v3 -> - begin - for l=0 to np-1 do - let v3_l = v3.(l) - and v1_l = v1.(l) - in - for k=0 to nq-1 do - v3_l.(k) <- v1_l.(k) +. v3_l.(k) - done - done; - Some v3 - end - with Exit -> p1 - in - try - let am = - try Powers.decr xyz angMom_a - with Invalid_argument _ -> raise Exit - in - if (axyz < 1) || (cxyz < 1) then - raise Exit; - let v = - vrr_v (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) - in - match (p2, v) with - | None, None -> None - | Some p2, None -> Some p2 - | _, Some v -> - begin - let p2 = - match p2 with - | None -> Array.make_matrix np nq 0. - | Some p2 -> p2 + begin + for l=0 to np-1 do + let v3_l = v3.(l) + and v1_l = v1.(l) + and p1_l = p1.(l) in - for l=0 to np-1 do - let fa = (float_of_int axyz) *. expo_inv_p.(l) *. 0.5 in - let p2_l = p2.(l) - and v_l = v.(l) - in - for k=0 to nq-1 do - p2_l.(k) <- p2_l.(k) -. fa *. expo_inv_q.(k) *. v_l.(k) - done - done; - Some p2 - end - with Exit -> p2 + for k=0 to nq-1 do + v3_l.(k) <- p1_l.(k) +. v1_l.(k) +. v3_l.(k) + done + done; + Some v3 + end + | Some p1, Some v1, None -> + begin + for l=0 to np-1 do + let v1_l = v1.(l) + and p1_l = p1.(l) + in + for k=0 to nq-1 do + p1_l.(k) <- v1_l.(k) +. p1_l.(k) + done + done; + Some p1 + end + | Some p1, None, Some v3 -> + begin + for l=0 to np-1 do + let v3_l = v3.(l) + and p1_l = p1.(l) + in + for k=0 to nq-1 do + p1_l.(k) <- p1_l.(k) +. v3_l.(k) + done + done; + Some p1 + end + | None , Some v1, Some v3 -> + begin + for l=0 to np-1 do + let v3_l = v3.(l) + and v1_l = v1.(l) + in + for k=0 to nq-1 do + v3_l.(k) <- v1_l.(k) +. v3_l.(k) + done + done; + Some v3 + end + in + if (axyz < 1) || (cxyz < 1) then p2 else + let am = Powers.decr xyz angMom_a in + let v = + vrr_v (m+1) am cm + in + match (p2, v) with + | None, None -> None + | Some p2, None -> Some p2 + | _, Some v -> + begin + let p2 = + match p2 with + | None -> Array.make_matrix np nq 0. + | Some p2 -> p2 + in + for l=0 to np-1 do + let fa = (float_of_int axyz) *. expo_inv_p.(l) *. 0.5 in + let p2_l = p2.(l) + and v_l = v.(l) + in + for k=0 to nq-1 do + p2_l.(k) <- p2_l.(k) -. fa *. expo_inv_q.(k) *. v_l.(k) + done + done; + Some p2 + end end in Zmap.add map_2d.(m) key result; result @@ -405,18 +369,17 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (** Horizontal recurrence relations *) - and hrr0_v angMom_a angMom_b angMom_c - totAngMom_a totAngMom_b totAngMom_c = + and hrr0_v angMom_a angMom_b angMom_c = - match totAngMom_b with + match angMom_b.tot with | 0 -> begin - match (totAngMom_a, totAngMom_c) with + match (angMom_a.tot, angMom_c.tot) with | (0,0) -> Array.fold_left (fun accu c -> - accu +. Array.fold_left (+.) 0. c) 0. zero_m_array.(0) + accu +. Array.fold_left (+.) 0. c) 0. zero_m_array.(0) | (_,_) -> begin - match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with + match vrr_v 0 angMom_a angMom_c with | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. end @@ -426,69 +389,60 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let ap = Powers.incr xyz angMom_a in let f = Coordinate.get xyz center_ab in let v1 = - match vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c with + match vrr_v 0 ap angMom_c with | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. in if (abs_float f < cutoff) then v1 else let v2 = - 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. + match vrr_v 0 angMom_a angMom_c with + | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix + | None -> 0. in v1 +. v2 *. f | _ -> let xyz = get_xyz angMom_b in let bxyz = Powers.get xyz angMom_b in - if (bxyz < 1) then 0. else + if (bxyz < 0) then 0. else let ap = Powers.incr xyz angMom_a in - let bm = - try Powers.decr xyz angMom_b - with Invalid_argument _ -> raise Exit - in - + let bm = Powers.decr xyz angMom_b in let h1 = - hrr0_v ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c + hrr0_v ap bm angMom_c in let f = Coordinate.get xyz center_ab in if abs_float f < cutoff then h1 else let h2 = - hrr0_v angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c + hrr0_v angMom_a bm angMom_c in h1 +. h2 *. f - and hrr_v angMom_a angMom_b angMom_c angMom_d - totAngMom_a totAngMom_b totAngMom_c totAngMom_d = + and hrr_v angMom_a angMom_b angMom_c angMom_d = - match (totAngMom_b, totAngMom_d) with - | (_,0) -> if totAngMom_b = 0 then + match (angMom_b.tot, angMom_d.tot) with + | (_,0) -> if angMom_b.tot = 0 then begin - match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with + match vrr_v 0 angMom_a angMom_c with | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix | None -> 0. end else - hrr0_v angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c + hrr0_v angMom_a angMom_b angMom_c | (_,_) -> let xyz = get_xyz angMom_d in let cp = Powers.incr xyz angMom_c in - let dm = - try Powers.decr xyz angMom_d - with Invalid_argument _ -> raise Exit - in + let dm = Powers.decr xyz angMom_d in let h1 = - hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) + hrr_v angMom_a angMom_b cp dm in let f = Coordinate.get xyz center_cd in 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) + hrr_v angMom_a angMom_b angMom_c dm in h1 +. f *. h2 in hrr_v angMom_a angMom_b angMom_c angMom_d - totAngMom_a totAngMom_b totAngMom_c totAngMom_d @@ -711,32 +665,27 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in (* Compute the integral class from the primitive shell quartet *) Array.iteri (fun i key -> - let (angMomA,angMomB,angMomC,angMomD) = - match Zkey.to_int_tuple ~kind:Zkey.Kind_12 key with + let (angMom_a,angMom_b,angMom_c,angMom_d) = + match Zkey.to_powers ~kind:Zkey.Kind_12 key with | Zkey.Twelve x -> x | _ -> assert false in try - if monocentric then - begin - 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; + if monocentric then + begin + if ( ((1 land angMom_a.x + angMom_b.x + angMom_c.x + angMom_d.x) =1) || + ((1 land angMom_a.y + angMom_b.y + angMom_c.y + angMom_d.y) =1) || + ((1 land angMom_a.z + angMom_b.z + angMom_c.z + angMom_d.z) =1) + ) then + raise NullQuartet + end; (* Schwartz screening *) if (np+nq> 24) then ( let schwartz_p = - let key = Zkey.of_int_tuple (Zkey.Twelve - (angMomA, angMomB, angMomA, angMomB) ) + let key = Zkey.of_powers (Zkey.Twelve + (angMom_a, angMom_b, angMom_a, angMom_b) ) in match schwartz_p with | None -> 1. @@ -744,8 +693,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in if schwartz_p < cutoff then raise NullQuartet; let schwartz_q = - let key = Zkey.of_int_tuple (Zkey.Twelve - (angMomC, angMomD, angMomC, angMomD) ) + let key = Zkey.of_powers (Zkey.Twelve + (angMom_c, angMom_d, angMom_c, angMom_d) ) in match schwartz_q with | None -> 1. @@ -755,10 +704,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q ); let integral = - hvrr_two_e_vector (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) + hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) + zero_m_array (expo_b, expo_d) (expo_inv_p, expo_inv_q) (shell_p.ContractedShellPair.center_ab, diff --git a/Utils/Angular_momentum.ml b/Utils/Angular_momentum.ml index 5f17a70..d7877a2 100644 --- a/Utils/Angular_momentum.ml +++ b/Utils/Angular_momentum.ml @@ -1,3 +1,5 @@ +open Powers + exception AngularMomentumError of string type t = @@ -69,37 +71,37 @@ let n_functions a = (** Returns an array of Zkeys corresponding to all possible angular momenta *) let zkey_array a = let keys_1d l = - let create_z (x,y,_) = - (x,y,l-(x+y)) + let create_z { x ; y ; _ } = + Powers.of_int_tuple (x,y,l-(x+y)) in let rec create_y accu xyz = - let (x,y,z) = xyz in + let { x ; y ; z } = xyz in match y with | 0 -> (create_z xyz)::accu | i -> let ynew = y-1 in - create_y ( (create_z xyz)::accu) (x,ynew,z) + create_y ( (create_z xyz)::accu) (Powers.of_int_tuple (x,ynew,z)) in let rec create_x accu xyz = - let (x,y,z) = xyz in + let { x ; y ; z } = xyz in match x with | 0 -> (create_y [] xyz)@accu | i -> let xnew = x-1 in let ynew = l-xnew in - create_x ((create_y [] xyz)@accu) (xnew, ynew, z) + create_x ((create_y [] xyz)@accu) (Powers.of_int_tuple (xnew, ynew, z)) in - create_x [] (l,0,0) + create_x [] (Powers.of_int_tuple (l,0,0)) |> List.rev in begin match a with | Singlet l1 -> - List.map (fun x -> Zkey.of_int_tuple (Zkey.Three x)) (keys_1d @@ to_int l1) + List.map (fun x -> Zkey.of_powers (Zkey.Three x)) (keys_1d @@ to_int l1) | Doublet (l1, l2) -> List.map (fun a -> List.map (fun b -> - Zkey.of_int_tuple (Zkey.Six (a,b))) (keys_1d @@ to_int l2) + Zkey.of_powers (Zkey.Six (a,b))) (keys_1d @@ to_int l2) ) (keys_1d @@ to_int l1) |> List.concat @@ -108,7 +110,7 @@ let zkey_array a = List.map (fun a -> List.map (fun b -> List.map (fun c -> - Zkey.of_int_tuple (Zkey.Nine (a,b,c))) (keys_1d @@ to_int l3) + Zkey.of_powers (Zkey.Nine (a,b,c))) (keys_1d @@ to_int l3) ) (keys_1d @@ to_int l2) |> List.concat ) (keys_1d @@ to_int l1) @@ -120,7 +122,7 @@ let zkey_array a = List.map (fun b -> List.map (fun c -> List.map (fun d -> - Zkey.of_int_tuple (Zkey.Twelve (a,b,c,d))) (keys_1d @@ to_int l4) + Zkey.of_powers (Zkey.Twelve (a,b,c,d))) (keys_1d @@ to_int l4) ) (keys_1d @@ to_int l3) |> List.concat ) (keys_1d @@ to_int l2) diff --git a/Utils/Powers.ml b/Utils/Powers.ml index 3fd8998..9d726a3 100644 --- a/Utils/Powers.ml +++ b/Utils/Powers.ml @@ -18,22 +18,15 @@ let get coord t = let incr coord t = match coord with - | Coordinate.X -> { t with x = t.x+1 ; tot = t.tot+1 } - | Coordinate.Y -> { t with y = t.y+1 ; tot = t.tot+1 } - | Coordinate.Z -> { t with z = t.z+1 ; tot = t.tot+1 } + | Coordinate.X -> let r = t.x+1 in { t with x = r ; tot = t.tot+1 } + | Coordinate.Y -> let r = t.y+1 in { t with y = r ; tot = t.tot+1 } + | Coordinate.Z -> let r = t.z+1 in { t with z = r ; tot = t.tot+1 } let decr coord t = - (* - let test _ = () - *) - let test x = - if x < 1 then - invalid_arg "Angular_momentum.Powers.decr"; - in match coord with - | Coordinate.X -> (test t.x ; { t with x = t.x-1 ; tot = t.tot-1 }) - | Coordinate.Y -> (test t.y ; { t with y = t.y-1 ; tot = t.tot-1 }) - | Coordinate.Z -> (test t.z ; { t with z = t.z-1 ; tot = t.tot-1 }) + | Coordinate.X -> let r = t.x-1 in { t with x = r ; tot = t.tot-1 } + | Coordinate.Y -> let r = t.y-1 in { t with y = r ; tot = t.tot-1 } + | Coordinate.Z -> let r = t.z-1 in { t with z = r ; tot = t.tot-1 } diff --git a/Utils/Powers.mli b/Utils/Powers.mli index da45bde..d4b05d2 100644 --- a/Utils/Powers.mli +++ b/Utils/Powers.mli @@ -2,6 +2,6 @@ type t = private { x: int ; y : int ; z : int ; tot : int } val of_int_tuple : int * int * int -> t val to_int_tuple : t -> int * int * int val get : Coordinate.axis -> t -> int -val incr : Coordinate.axis -> t -> t +val incr : Coordinate.axis -> t -> t val decr : Coordinate.axis -> t -> t diff --git a/Utils/Zkey.ml b/Utils/Zkey.ml index b15bd56..f5946f5 100644 --- a/Utils/Zkey.ml +++ b/Utils/Zkey.ml @@ -1,12 +1,11 @@ +open Powers + (** Key for hastables that contain tuples of integers encoded in small integers *) type kind_array = | Kind_3 | Kind_6 | Kind_12 | Kind_9 -| Kind_4 -| Kind_2 -| Kind_1 type t = { @@ -40,32 +39,23 @@ let of_int_array ~kind a = | Kind_9 -> of_int a.(0) << a.(1) << a.(2) << a.(3) << a.(4) << a.(5) <| a.(6) << a.(7) << a.(8) - | Kind_4 -> of_int a.(0) <+ a.(1) <+ a.(2) <+ a.(3) - | Kind_2 -> of_int a.(0) <+ a.(1) - | Kind_1 -> of_int a.(0) type kind = - | One of (int) - | Two of (int*int) - | Three of (int*int*int) - | Four of ((int*int)*(int*int)) - | Six of ((int*int*int)*(int*int*int)) - | Nine of ((int*int*int)*(int*int*int)*(int*int*int)) - | Twelve of ((int*int*int)*(int*int*int)*(int*int*int)*(int*int*int)) + | Three of Powers.t + | Six of (Powers.t * Powers.t) + | Nine of (Powers.t * Powers.t * Powers.t) + | Twelve of (Powers.t * Powers.t * Powers.t * Powers.t) -let of_int_tuple a = +let of_powers a = match a with - | One (a) -> of_int a - | Two (a,b) -> of_int a <+ b - | Three (a,b,c) -> of_int a <+ b <+ c - | Four ((a,b),(c,d)) -> of_int a <+ b <+ c <+ d - | Six ((a,b,c),(d,e,f)) -> + | Three { x=a ; y=b ; z=c ; _ } -> of_int a <+ b <+ c + | Six ({ x=a ; y=b ; z=c },{ x=d ; y=e ; z=f }) -> of_int a << b << c << d << e << f - | Twelve ((a,b,c),(d,e,f),(g,h,i),(j,k,l)) -> + | Twelve ({ x=a ; y=b ; z=c },{ x=d ; y=e ; z=f },{ x=g ; y=h ; z=i },{ x=j ; y=k ; z=l }) -> of_int a << b << c << d << e << f <| g << h << i << j << k << l - | Nine ((a,b,c),(d,e,f),(g,h,i)) -> + | Nine ({ x=a ; y=b ; z=c },{ x=d ; y=e ; z=f },{ x=g ; y=h ; z=i }) -> of_int a << b << c << d << e << f <| g << h << i @@ -118,80 +108,59 @@ let to_int_array ~kind { left ; right } = mask10 land right |] - | Kind_4 -> [| - mask15 land (right lsr 45) ; - mask15 land (right lsr 30) ; - mask15 land (right lsr 15) ; - mask15 land right - |] - - | Kind_2 -> [| - mask15 land (right lsr 15) ; - mask15 land right - |] - - | Kind_1 -> [| right |] - (** Transform the Zkey into an int tuple *) -let to_int_tuple ~kind { left ; right } = +let to_powers ~kind { left ; right } = match kind with - | Kind_3 -> Three ( + | Kind_3 -> Three (Powers.of_int_tuple ( mask15 land (right lsr 30) , mask15 land (right lsr 15) , mask15 land right - ) + )) - | Kind_6 -> Six ( + | Kind_6 -> Six (Powers.of_int_tuple ( mask10 land (right lsr 50) , mask10 land (right lsr 40) , mask10 land (right lsr 30)), + Powers.of_int_tuple ( mask10 land (right lsr 20) , mask10 land (right lsr 10) , mask10 land right ) ) - | Kind_12 -> Twelve ( + | Kind_12 -> Twelve (Powers.of_int_tuple ( mask10 land (left lsr 50) , mask10 land (left lsr 40) , mask10 land (left lsr 30)), + Powers.of_int_tuple ( mask10 land (left lsr 20) , mask10 land (left lsr 10) , mask10 land left ) , + Powers.of_int_tuple ( mask10 land (right lsr 50) , mask10 land (right lsr 40) , mask10 land (right lsr 30)), + Powers.of_int_tuple ( mask10 land (right lsr 20) , mask10 land (right lsr 10) , mask10 land right ) ) - | Kind_9 -> Nine ( + | Kind_9 -> Nine (Powers.of_int_tuple ( mask10 land (left lsr 20) , mask10 land (left lsr 10) , mask10 land left ) , + Powers.of_int_tuple ( mask10 land (right lsr 50) , mask10 land (right lsr 40) , mask10 land (right lsr 30)), + Powers.of_int_tuple ( mask10 land (right lsr 20) , mask10 land (right lsr 10) , mask10 land right ) ) - | Kind_4 -> Four ( - ( mask15 land (right lsr 45) , - mask15 land (right lsr 30)), - ( mask15 land (right lsr 15) , - mask15 land right ) - ) - - | Kind_2 -> Two ( - mask15 land (right lsr 15) , - mask15 land right - ) - - | Kind_1 -> One right let hash = Hashtbl.hash @@ -220,15 +189,11 @@ let to_string ~kind { left ; right } = (* let debug () = - let k2 = of_int_array Kind_2 [| 1 ; 2 |] and k3 = of_int_array Kind_3 [| 1 ; 2 ; 3 |] - and k4 = of_int_array Kind_4 [| 1 ; 2 ; 3; 4 |] and k6 = of_int_array Kind_6 [| 1 ; 2 ; 3; 4 ; 5; 6|] and k12 = of_int_array Kind_12 [| 1 ; 2 ; 3; 4 ; 5; 6 ; 7 ; 8 ; 9 ; 10 ; 11; 12|] in - print_endline @@ to_string Kind_2 k2 ; print_endline @@ to_string Kind_3 k3 ; - print_endline @@ to_string Kind_4 k4 ; print_endline @@ to_string Kind_6 k6 ; print_endline @@ to_string Kind_12 k12