From 8993c168d6e1238d8f3b2e869b48d6ae9e2be8b4 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 2 Feb 2018 10:49:40 +0100 Subject: [PATCH] Introducing tuples --- Basis/TwoElectronRR.ml | 171 +++++++++++++++++++------------ Basis/TwoElectronRRVectorized.ml | 32 +++--- 2 files changed, 124 insertions(+), 79 deletions(-) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index e83a0ae..7ad24c5 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -39,41 +39,40 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (** Vertical recurrence relations *) let rec vrr0 angMom_a totAngMom_a = - (* 1_137_164 *) if debug then - Printf.printf "vrr0: %d : %d %d %d\n" totAngMom_a angMom_a.(0) angMom_a.(1) angMom_a.(2); + (* 1_137_164 *) + if debug then + begin + let (x,y,z) = angMom_a in + Printf.printf "vrr0: %d : %d %d %d\n" totAngMom_a x y z + end; match totAngMom_a with | 0 -> (* 66_288 *) zero_m_array | _ -> (* 1_070_876 *) let maxsze = maxm+1 in - let key = Zkey.of_int_tuple (Zkey.Three - (angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1) ) + let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in try Zmap.find map_1d key with | Not_found -> let result = - let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in - let xyz = + let am, amm, amxyz, xyz = match angMom_a with - | [|_;0;0|] -> (* 28_336 *) 0 - | [|_;_;0|] -> (* 52_221 *) 1 - | _ -> (* 87_215 *) 2 + | (x,0,0) -> (* 28_336 *) (x-1,0,0),(x-2,0,0), x-1, 0 + | (x,y,0) -> (* 52_221 *) (x,y-1,0),(x,y-2,0), y-1, 1 + | (x,y,z) -> (* 87_215 *) (x,y,z-1),(x,y,z-2), z-1, 2 in - am.(xyz) <- am.(xyz) - 1; - if am.(xyz) < 0 then empty else + if amxyz < 0 then empty else let v1 = vrr0 am (totAngMom_a-1) in let f1 = expo_inv_p *. (Coordinate.coord center_pq xyz) and f2 = expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz) in - if am.(xyz) < 1 then + if amxyz < 1 then Array.init maxsze (fun m -> (* 544_860 *) if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) ) else - let f3 = (float_of_int am.(xyz)) *. expo_inv_p *. 0.5 in - let amm = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in - let () = amm.(xyz) <- amm.(xyz) - 2 in + let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in let v3 = vrr0 amm (totAngMom_a-2) in @@ -89,32 +88,46 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and vrr angMom_a angMom_c totAngMom_a totAngMom_c = - (* 11_580_843 *) if debug then - Printf.printf "vrr : %d %d : %d %d %d %d %d %d\n" totAngMom_a totAngMom_c angMom_a.(0) angMom_a.(1) angMom_a.(2) angMom_c.(0) angMom_c.(1) angMom_c.(2); + (* 11_580_843 *) + 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 + angMom_ax angMom_ay angMom_az angMom_cx angMom_cy angMom_cz + end; match (totAngMom_a, totAngMom_c) with | (i,0) -> (* 959_629 *) if (i>0) then vrr0 angMom_a totAngMom_a else zero_m_array | (_,_) -> (* 10_621_214 *) let maxsze = maxm+1 in - let key = Zkey.of_int_tuple (Zkey.Six - ((angMom_a.(0)+1, angMom_a.(1)+1, angMom_a.(2)+1), - (angMom_c.(0)+1, angMom_c.(1)+1, angMom_c.(2)+1)) ) - - in + let key = Zkey.of_int_tuple (Zkey.Six (angMom_a, angMom_c) ) in try Zmap.find map_2d key with | Not_found -> let result = - let cm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] in - let xyz = + let am, cm, cmm, axyz, cxyz, 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|] -> (* 321_984 *) 0 - | [|_;_;0|] -> (* 612_002 *) 1 - | _ -> (* 1_067_324 *) 2 + | (_,0,0) -> (* 321_984 *) + (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, 0 + | (_,_,0) -> (* 612_002 *) + (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 + | _ -> (* 1_067_324 *) + (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, 2 in - cm.(xyz) <- cm.(xyz) - 1; - if cm.(xyz) < 0 then empty else + if cxyz < 1 then empty else let f1 = -. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz) in @@ -130,14 +143,12 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (* 9_232_029 *) f1 *. v1.(m) -. (if m = maxm then 0. else f2 *. v1.(m+1)) ) in let result = - if cm.(xyz) < 1 then result else + if cxyz < 2 then result else let f3 = - (float_of_int cm.(xyz)) *. expo_inv_q *. 0.5 + (float_of_int (cxyz-1)) *. expo_inv_q *. 0.5 in if (abs_float f3 < cutoff) && (abs_float (f3 *. expo_inv_q) < cutoff) then result else ( - let cmm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] in - cmm.(xyz) <- cmm.(xyz) - 2; let v3 = vrr angMom_a cmm totAngMom_a (totAngMom_c-2) in @@ -146,11 +157,9 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) ) in let result = - let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] in - am.(xyz) <- am.(xyz) - 1; - if am.(xyz) lor cm.(xyz) < 0 then result else + if axyz < 1 || cxyz < 1 then result else let f5 = - (float_of_int angMom_a.(xyz)) *. expo_inv_p *. expo_inv_q *. 0.5 + (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 in if (abs_float f5 < cutoff) then result else let v5 = @@ -171,15 +180,30 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c = - (* 8_448_486 *) if debug then - Printf.printf "hrr0: %d %d %d : %d %d %d %d %d %d %d %d %d\n" totAngMom_a totAngMom_b totAngMom_c 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); + (* 8_448_486 *) + if debug then + begin + let angMom_ax, angMom_ay, angMom_az = angMom_a + 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 -> (* 0 *) (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0) | 1 -> - (* 5_045_008 *) 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 angMom_ax, angMom_ay, angMom_az = angMom_a in + (* 5_045_008 *) + let ap, xyz = + match angMom_b with + | (1,_,_) -> (angMom_ax+1,angMom_ay,angMom_az), 0 + | (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), 1 + | (_,_,_) -> (angMom_ax,angMom_ay,angMom_az+1), 2 + in let v1 = vrr ap angMom_c (totAngMom_a+1) totAngMom_c in @@ -192,17 +216,22 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) in v1.(0) +. f2 *. v2.(0) | _ -> - (* 3_403_478 *) let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] - and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |] - and xyz = - match angMom_b with - | [|_;0;0|] -> (* 677_315 *) 0 - | [|_;_;0|] -> (* 1_136_646 *) 1 - | _ -> (* 1_589_517 *) 2 + let angMom_ax, angMom_ay, angMom_az = angMom_a + and angMom_bx, angMom_by, angMom_bz = angMom_b in + (* 3_403_478 *) + let bxyz, xyz = + match angMom_b with + | (1,_,_) -> angMom_bx, 0 + | (_,1,_) -> angMom_by, 1 + | (_,_,_) -> angMom_bz, 2 in - ap.(xyz) <- ap.(xyz) + 1; - bm.(xyz) <- bm.(xyz) - 1; - if (bm.(xyz) < 0) then 0. else + if (bxyz < 1) then 0. else + let ap, bm = + match xyz with + | 0 -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz) + | 1 -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz) + | _ -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1) + in let h1 = hrr0 ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c in @@ -218,8 +247,20 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and hrr angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d = - (* 7_738_602 *) if debug then - Printf.printf "hrr : %d %d %d %d : %d %d %d %d %d %d %d %d %d %d %d %d\n" totAngMom_a totAngMom_b totAngMom_c totAngMom_d 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) angMom_d.(0) angMom_d.(1) angMom_d.(2); + (* 7_738_602 *) + if debug then + begin + let angMom_ax, angMom_ay, angMom_az = angMom_a in + let angMom_bx, angMom_by, angMom_bz = angMom_b in + let angMom_cx, angMom_cy, angMom_cz = angMom_c in + let angMom_dx, angMom_dy, angMom_dz = angMom_d in + Printf.printf "hrr : %d %d %d %d : %d %d %d %d %d %d %d %d %d %d %d %d\n" + totAngMom_a totAngMom_b totAngMom_c totAngMom_d + angMom_ax angMom_ay angMom_az + angMom_bx angMom_by angMom_bz + angMom_cx angMom_cy angMom_cz + angMom_dx angMom_dy angMom_dz + end; match (totAngMom_b, totAngMom_d) with | (_,0) -> (* 3_608_781 *) if (totAngMom_b = 0) then @@ -227,16 +268,15 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) else hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c | (_,_) -> - (* 4_130_325 *) let cp = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] - and dm = [| angMom_d.(0) ; angMom_d.(1) ; angMom_d.(2) |] - and xyz = + (* 4_130_325 *) + 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|] -> (* 1_524_451 *) 0 - | [|_;_;0|] -> (* 1_302_937 *) 1 - | _ -> (* 1_302_937 *) 2 + | (_,0,0) -> (* 1_524_451 *) (angMom_cx+1, angMom_cy, angMom_cz), (angMom_dx-1, angMom_dy, angMom_dz), 0 + | (_,_,0) -> (* 1_302_937 *) (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), 1 + | _ -> (* 1_302_937 *) (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), 2 in - cp.(xyz) <- cp.(xyz) + 1; - dm.(xyz) <- dm.(xyz) - 1; let h1 = hrr angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) in @@ -248,7 +288,12 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) h1 +. f2 *. h2 in - hrr angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d + hrr + (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)) + (angMom_d.(0),angMom_d.(1),angMom_d.(2)) + totAngMom_a totAngMom_b totAngMom_c totAngMom_d diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index f17ecae..307370e 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -53,10 +53,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] and amm = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] and xyz = - match angMom_a with - | [|0;0;_|] -> 2 - | [|0;_;_|] -> 1 - | _ -> 0 + match (angMom_a.(0),angMom_a.(1),angMom_a.(2)) with + | (0,0,_) -> 2 + | (0,_,_) -> 1 + | _ -> 0 in am.(xyz) <- am.(xyz) - 1; amm.(xyz) <- amm.(xyz) - 2; @@ -112,10 +112,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) and cm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] and cmm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] and xyz = - match angMom_c with - | [|0;0;_|] -> 2 - | [|0;_;_|] -> 1 - | _ -> 0 + match (angMom_c.(0),angMom_c.(1),angMom_c.(2)) with + | (0,0,_) -> 2 + | (0,_,_) -> 1 + | _ -> 0 in am.(xyz) <- am.(xyz) - 1; cm.(xyz) <- cm.(xyz) - 1; @@ -217,10 +217,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |] and xyz = - match angMom_b with - | [|0;0;_|] -> 2 - | [|0;_;_|] -> 1 - | _ -> 0 + match (angMom_b.(0),angMom_b.(1),angMom_b.(2)) with + | (0,0,_) -> 2 + | (0,_,_) -> 1 + | _ -> 0 in ap.(xyz) <- ap.(xyz) + 1; bm.(xyz) <- bm.(xyz) - 1; @@ -244,10 +244,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let cp = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] and dm = [| angMom_d.(0) ; angMom_d.(1) ; angMom_d.(2) |] and xyz = - match angMom_d with - | [|0;0;_|] -> 2 - | [|0;_;_|] -> 1 - | _ -> 0 + match (angMom_d.(0),angMom_d.(1),angMom_d.(2)) with + | (0,0,_) -> 2 + | (0,_,_) -> 1 + | _ -> 0 in cp.(xyz) <- cp.(xyz) + 1; dm.(xyz) <- dm.(xyz) - 1;