From 7bcc910edf665213cb570af0fb496c5d255bc6e2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 9 Feb 2018 10:49:27 +0100 Subject: [PATCH] Swapping of 2 electrons for HRR --- Basis/TwoElectronRR.ml | 491 +++++++++++++++++++++-------------------- Utils/Coordinate.ml | 2 + Utils/Coordinate.mli | 1 + 3 files changed, 255 insertions(+), 239 deletions(-) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 9bc6986..2a24c48 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -8,13 +8,13 @@ let cutoff2 = cutoff *. cutoff exception NullQuartet (** Horizontal and Vertical Recurrence Relations (HVRR) *) -let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) +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) (expo_b, expo_d) (expo_inv_p, expo_inv_q) (center_ab, center_cd, center_pq) - map_1d map_2d + map_1d map_2d map_1d' map_2d' = let maxsze = maxm+1 in @@ -23,276 +23,287 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and totAngMom_c = Angular_momentum.to_int totAngMom_c_in and totAngMom_d = Angular_momentum.to_int totAngMom_d_in in - let maxm = totAngMom_a + totAngMom_b + totAngMom_c + totAngMom_d in - let empty = Array.make (maxm+1) 0. - 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; - let (x,y,z) = angMom_b in Printf.printf "%d %d %d\n" x y z; - let (x,y,z) = angMom_c in Printf.printf "%d %d %d\n" x y z; - let (x,y,z) = angMom_d in Printf.printf "%d %d %d\n" x y z; - Printf.printf "%f %f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n" expo_b expo_d - expo_inv_p expo_inv_q - (Coordinate.coord center_ab 0) (Coordinate.coord center_ab 1) (Coordinate.coord center_ab 2) - (Coordinate.coord center_cd 0) (Coordinate.coord center_cd 1) (Coordinate.coord center_cd 2) - (Coordinate.coord center_pq 0) (Coordinate.coord center_pq 1) (Coordinate.coord center_pq 2) - end; - (** Vertical recurrence relations *) - let rec vrr0 angMom_a totAngMom_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 - end; + (* Swap electrons 1 and 2 so that the max angular momentum is on 1 *) + if (totAngMom_a+totAngMom_b < totAngMom_c+totAngMom_d) 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) + (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. + 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; + let (x,y,z) = angMom_b in Printf.printf "%d %d %d\n" x y z; + let (x,y,z) = angMom_c in Printf.printf "%d %d %d\n" x y z; + let (x,y,z) = angMom_d in Printf.printf "%d %d %d\n" x y z; + Printf.printf "%f %f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n" expo_b expo_d + expo_inv_p expo_inv_q + (Coordinate.coord center_ab 0) (Coordinate.coord center_ab 1) (Coordinate.coord center_ab 2) + (Coordinate.coord center_cd 0) (Coordinate.coord center_cd 1) (Coordinate.coord center_cd 2) + (Coordinate.coord center_pq 0) (Coordinate.coord center_pq 1) (Coordinate.coord center_pq 2) + end; - match totAngMom_a with - | 0 -> zero_m_array - | _ -> - let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in + (** Vertical recurrence relations *) + let rec vrr0 angMom_a totAngMom_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 + end; - 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, 0 - | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1 - | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2 - in - if amxyz < 0 then empty else - let f1 = expo_inv_p *. (Coordinate.coord center_pq xyz) - and f2 = expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz) + match totAngMom_a with + | 0 -> zero_m_array + | _ -> + let key = Zkey.of_int_tuple (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, 0 + | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1 + | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2 in - if amxyz < 1 then - let v1 = - vrr0 am (totAngMom_a-1) + if amxyz < 0 then empty else + let f1 = expo_inv_p *. (Coordinate.coord center_pq xyz) + and f2 = expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz) in - Array.init maxsze (fun m -> - if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) ) - else - let v3 = - vrr0 amm (totAngMom_a-2) - in - let v1 = - vrr0 am (totAngMom_a-1) - in - let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in - Array.init maxsze (fun m -> - (if m = maxm then 0. else - (f1 *. v1.(m+1) ) -. f2 *. v1.(m) ) - +. f3 *. (v3.(m) +. if m = maxm then 0. else - expo_inv_p *. v3.(m+1)) - ) - in Zmap.add map_1d key result; - result + if amxyz < 1 then + let v1 = + vrr0 am (totAngMom_a-1) + in + Array.init maxsze (fun m -> + if m = maxm then 0. else (f1 *. v1.(m+1) ) -. f2 *. v1.(m) ) + else + let v3 = + vrr0 amm (totAngMom_a-2) + in + let v1 = + vrr0 am (totAngMom_a-1) + in + let f3 = (float_of_int amxyz) *. expo_inv_p *. 0.5 in + Array.init maxsze (fun m -> + (if m = maxm then 0. else + (f1 *. v1.(m+1) ) -. f2 *. v1.(m) ) + +. f3 *. (v3.(m) +. if m = maxm then 0. else + expo_inv_p *. v3.(m+1)) + ) + in Zmap.add map_1d key result; + result - and vrr angMom_a angMom_c totAngMom_a totAngMom_c = + and vrr angMom_a angMom_c totAngMom_a totAngMom_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 - angMom_ax angMom_ay angMom_az angMom_cx angMom_cy angMom_cz - end; + 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) -> if (i>0) then - vrr0 angMom_a totAngMom_a + match (totAngMom_a, totAngMom_c) with + | (i,0) -> if (i>0) then + vrr0 angMom_a totAngMom_a (* OneElectronRR.hvrr_one_e (angMom_a, angMom_b) (totAngMom_a_in, totAngMom_b_in) (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 + else zero_m_array + | (_,_) -> + 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 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) -> (* 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-1, 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 - | _ -> (* 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-1, 2 - in - if cmxyz < 0 then empty else - let f1 = - -. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz) + 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) -> (* 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-1, 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 + | _ -> (* 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-1, 2 in - let f2 = - expo_inv_q *. (Coordinate.coord center_pq xyz) - in - let result = empty in - let result = - if ( (abs_float f1 < cutoff) && (abs_float f2 < cutoff) ) then result else - let v1 = - vrr angMom_a cm totAngMom_a (totAngMom_c-1) - in - Array.init maxsze (fun m -> result.(m) +. - f1 *. v1.(m) -. (if m = maxm then 0. else f2 *. v1.(m+1)) ) - in - let result = - if cmxyz < 1 then result else - let f3 = - (float_of_int cmxyz) *. expo_inv_q *. 0.5 - in - if (abs_float f3 < cutoff) && (abs_float (f3 *. expo_inv_q) < cutoff) then result else - ( - let v3 = - vrr angMom_a cmm totAngMom_a (totAngMom_c-2) - in - Array.init maxsze (fun m -> result.(m) +. - f3 *. (v3.(m) +. (if m=maxm then 0. else expo_inv_q *. v3.(m+1)) )) - ) - in - let result = - if (axyz < 1) || (cmxyz < 0) then result else - let f5 = - (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 - in - if (abs_float f5 < cutoff) then result else - let v5 = - vrr am cm (totAngMom_a-1) (totAngMom_c-1) + if cmxyz < 0 then empty else + let f1 = + -. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz) + in + let f2 = + expo_inv_q *. (Coordinate.coord center_pq xyz) + in + let result = empty in + let result = + if ( (abs_float f1 < cutoff) && (abs_float f2 < cutoff) ) then result else + let v1 = + vrr angMom_a cm totAngMom_a (totAngMom_c-1) in - Array.init (maxsze) (fun m -> - result.(m) -. (if m = maxm then 0. else f5 *. v5.(m+1))) - in - result - in Zmap.add map_2d key result; - result + Array.init maxsze (fun m -> result.(m) +. + f1 *. v1.(m) -. (if m = maxm then 0. else f2 *. v1.(m+1)) ) + in + let result = + if cmxyz < 1 then result else + let f3 = + (float_of_int cmxyz) *. expo_inv_q *. 0.5 + in + if (abs_float f3 < cutoff) && (abs_float (f3 *. expo_inv_q) < cutoff) then result else + ( + let v3 = + vrr angMom_a cmm totAngMom_a (totAngMom_c-2) + in + Array.init maxsze (fun m -> result.(m) +. + f3 *. (v3.(m) +. (if m=maxm then 0. else expo_inv_q *. v3.(m+1)) )) + ) + in + let result = + if (axyz < 1) || (cmxyz < 0) then result else + let f5 = + (float_of_int axyz) *. expo_inv_p *. expo_inv_q *. 0.5 + in + if (abs_float f5 < cutoff) then result else + let v5 = + vrr am cm (totAngMom_a-1) (totAngMom_c-1) + in + Array.init (maxsze) (fun m -> + result.(m) -. (if m = maxm then 0. else f5 *. v5.(m+1))) + in + result + in Zmap.add map_2d key result; + result - (** Horizontal recurrence relations *) - and hrr0 angMom_a angMom_b angMom_c - totAngMom_a totAngMom_b totAngMom_c = + (** Horizontal recurrence relations *) + and hrr0 angMom_a angMom_b angMom_c + totAngMom_a totAngMom_b totAngMom_c = - 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; + 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 -> (vrr angMom_a angMom_c totAngMom_a totAngMom_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), 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 - let f2 = - (Coordinate.coord center_ab xyz) - in - if (abs_float f2 < cutoff) then v1.(0) else - let v2 = - vrr angMom_a angMom_c totAngMom_a totAngMom_c + match totAngMom_b with + | 0 -> (vrr angMom_a angMom_c totAngMom_a totAngMom_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), 0 + | (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), 1 + | _ -> (angMom_ax,angMom_ay,angMom_az+1), 2 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, 0 - | (_,_,0) -> angMom_by, 1 - | (_,_,_) -> angMom_bz, 2 - in - 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 + let v1 = + vrr ap angMom_c (totAngMom_a+1) totAngMom_c in let f2 = (Coordinate.coord center_ab xyz) in + if (abs_float f2 < cutoff) then v1.(0) else + let v2 = + vrr angMom_a angMom_c totAngMom_a totAngMom_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, 0 + | (_,_,0) -> angMom_by, 1 + | (_,_,_) -> angMom_bz, 2 + in + 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 + let f2 = + (Coordinate.coord center_ab xyz) + in + if (abs_float f2 < cutoff) then h1 else + let h2 = + hrr0 angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c + in + h1 +. f2 *. h2 + + and hrr angMom_a angMom_b angMom_c angMom_d + totAngMom_a totAngMom_b totAngMom_c totAngMom_d = + + 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) -> if (totAngMom_b = 0) then + (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0) + else + hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_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), 0 + | (_,_,0) -> (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), 1 + | _ -> (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), 2 + in + let h1 = + hrr angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) + in + let f2 = Coordinate.coord center_cd xyz in if (abs_float f2 < cutoff) then h1 else let h2 = - hrr0 angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c + hrr angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) in h1 +. f2 *. h2 - and hrr angMom_a angMom_b angMom_c angMom_d - totAngMom_a totAngMom_b totAngMom_c totAngMom_d = - - 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) -> if (totAngMom_b = 0) then - (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0) - else - hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_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), 0 - | (_,_,0) -> (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), 1 - | _ -> (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), 2 - in - let h1 = - hrr angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) - in - let f2 = Coordinate.coord center_cd xyz 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) - in - h1 +. f2 *. h2 - - in - hrr angMom_a angMom_b angMom_c angMom_d - totAngMom_a totAngMom_b totAngMom_c totAngMom_d + in + hrr angMom_a angMom_b angMom_c angMom_d + totAngMom_a totAngMom_b totAngMom_c totAngMom_d @@ -362,6 +373,8 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let d = shell_q.(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_q = shell_q.(cd).ShellPair.norm_coef_scale in let norm_coef_scale = Array.map (fun v1 -> @@ -435,7 +448,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d) (shell_p.(ab).ShellPair.expo_inv, shell_q.(cd).ShellPair.expo_inv) (shell_p.(ab).ShellPair.center_ab, shell_q.(cd).ShellPair.center_ab, center_pq) - map_1d map_2d + map_1d map_2d map_1d' map_2d' in contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral with NullQuartet -> () diff --git a/Utils/Coordinate.ml b/Utils/Coordinate.ml index 797ee0e..eaabfa5 100644 --- a/Utils/Coordinate.ml +++ b/Utils/Coordinate.ml @@ -54,6 +54,8 @@ let (|-), (|+) = ) +let neg a = -1. |. a + let rec dot p q = match (p,q) with | Bohr (x,y,z), Bohr (x',y',z') -> x*.x' +. y*.y' +. z*.z' diff --git a/Utils/Coordinate.mli b/Utils/Coordinate.mli index 0e1599c..13df088 100644 --- a/Utils/Coordinate.mli +++ b/Utils/Coordinate.mli @@ -7,6 +7,7 @@ val of_3_floats : float -> float -> float -> [< `Angstrom | `Bohr ] -> t val ( |. ) : float -> t -> t val ( |- ) : t -> t -> t val ( |+ ) : t -> t -> t +val neg : t -> t val dot : t -> t -> float val norm : t -> float val to_string : t -> string