diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 23a2fb1..58a99de 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -196,31 +196,16 @@ let to_file ~filename basis = Printf.printf "Computed %d non-zero ERIs in %f seconds\n" !inn (Unix.gettimeofday () -. t0); (* Print ERIs *) - for i=0 to (Array.length basis) - 1 do - print_int basis.(i).Contracted_shell.indice ; print_newline (); - for j=0 to i do - for k=0 to i do - for l=0 to k do - (* Write the data in the output file *) - Array.iteri (fun i_c _ -> - let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in - Array.iteri (fun j_c _ -> - let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in - Array.iteri (fun k_c _ -> - let k_c = basis.(k).Contracted_shell.indice + k_c + 1 in - Array.iteri (fun l_c _ -> - let l_c = basis.(l).Contracted_shell.indice + l_c + 1 in - let value = eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} in - if (value <> 0.) then - Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" - i_c k_c j_c l_c value; - ) basis.(l).Contracted_shell.powers - ) basis.(k).Contracted_shell.powers - ) basis.(j).Contracted_shell.powers - ) basis.(i).Contracted_shell.powers; - done; - done; + for i_c=1 to (Genarray.nth_dim eri_array 0) do + for j_c=1 to (Genarray.nth_dim eri_array 2) do + for k_c=1 to (Genarray.nth_dim eri_array 1) do + for l_c=1 to (Genarray.nth_dim eri_array 3) do + let value = eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} in + if (value <> 0.) then + Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i_c k_c j_c l_c value; + done; done; + done; done; Printf.printf "In: %d Out:%d\n" !inn !out ; close_out oc diff --git a/Basis/Shell_pair.ml b/Basis/Shell_pair.ml index 0807f9c..2f2ae58 100644 --- a/Basis/Shell_pair.ml +++ b/Basis/Shell_pair.ml @@ -15,6 +15,7 @@ type t = { j : int; shell_a : Contracted_shell.t; shell_b : Contracted_shell.t; + monocentric : bool } exception Null_contribution @@ -86,7 +87,10 @@ let create_array ?cutoff p_a p_b = let center_a = Coordinate.(center |- Contracted_shell.center p_a) in - Some { i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; norm_coef_scale } + let monocentric = + Contracted_shell.center p_a = Contracted_shell.center p_b + in + Some { i ; j ; shell_a=p_a ; shell_b=p_b ; norm_coef ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq ; norm_coef_scale ; monocentric } with | Null_contribution -> None ) diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 0665480..6eff461 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -31,11 +31,11 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) Printf.printf "%d %d %d\n" angMom_c.(0) angMom_c.(1) angMom_c.(2) ; Printf.printf "%d %d %d\n" angMom_d.(0) angMom_d.(1) angMom_d.(2) ; Printf.printf "%f %f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n" expo_b expo_d - end 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); + (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 = @@ -53,15 +53,12 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) try Zmap.find map_1d key with | Not_found -> let result = - let am, amm, amxyz, xyz = - let x, y, z = angMom_a in - if (z > 0) then - (x,y,z-1),(x,y,z-2), z-1, 2 - else if (y > 0) then - (x,y-1,0),(x,y-2,0), y-1, 1 - else - (x-1,0,0),(x-2,0,0), x-1, 0 - in + 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 v1 = vrr0 am (totAngMom_a-1) @@ -125,7 +122,9 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (angMom_cx-2, angMom_cy, angMom_cz), angMom_ax,angMom_cx, 0 in + (* if cxyz < 1 then empty else + *) let f1 = -. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz) in @@ -364,6 +363,12 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q |> Array.to_list |> Array.concat in + (* + let monocentric = + shell_p.(ab).Shell_pair.monocentric && + shell_q.(cd).Shell_pair.monocentric + in + *) (* Compute the integral class from the primitive shell quartet *) class_indices |> Array.iteri (fun i key -> @@ -375,6 +380,16 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q [| a.(9) ; a.(10) ; a.(11) |] ) in try + (* + if monocentric then + begin + if ( ((1 land (a.(0) + a.(3) + a.(6) + a.( 9)))=1) || + ((1 land (a.(1) + a.(4) + a.(7) + a.(10)))=1) || + ((1 land (a.(2) + a.(5) + a.(8) + a.(11)))=1) + ) then + raise NullQuartet + end; + *) (* Schwartz screening *) (* let schwartz_p = @@ -408,7 +423,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 = - (* 2_011_273 *) hvrr_two_e (angMomA, angMomB, angMomC, angMomD) + 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) diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index e7965fe..5efdf0b 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -103,25 +103,27 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result = begin 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 - if (angMom_cz > 0) then - (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 - else if (angMom_cy > 0) then - (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 + let (aax, aay, aaz) = angMom_a + and (acx, acy, acz) = angMom_c in + if (acz > 0) then + (aax, aay, aaz-1), + (acx, acy, acz-1), + (acx, acy, acz-2), + aaz, acz, 2 + else if (acy > 0) then + (aax, aay-1,aaz), + (acx, acy-1,acz), + (acx, acy-2,acz), + aay,acy, 1 else - (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 + (aax-1,aay,aaz), + (acx-1,acy,acz), + (acx-2,acy,acz), + aax,acx, 0 in + (* if cxyz < 1 then empty else + *) let f1 = Array.init ncoef (fun k -> expo_d.(k) *. expo_inv_q.(k) *. @@ -199,12 +201,12 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) | (_,_) -> vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c end | 1 -> - let angMom_ax, angMom_ay, angMom_az = angMom_a in + let (aax, aay, aaz) = angMom_a in let ap, xyz = match angMom_b with - | (_,_,1) -> (angMom_ax,angMom_ay,angMom_az+1), 2 - | (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), 1 - | (_,_,_) -> (angMom_ax+1,angMom_ay,angMom_az), 0 + | (_,_,1) -> (aax,aay,aaz+1), 2 + | (_,1,_) -> (aax,aay+1,aaz), 1 + | (_,_,_) -> (aax+1,aay,aaz), 0 in let f = Coordinate.coord center_ab xyz in let v1 = @@ -216,20 +218,20 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) in Array.map2 (fun v1 v2 -> v1 +. v2 *. f) v1 v2 | _ -> - let angMom_ax, angMom_ay, angMom_az = angMom_a - and angMom_bx, angMom_by, angMom_bz = angMom_b in + let (aax, aay, aaz) = angMom_a + and (abx, aby, abz) = angMom_b in let bxyz, xyz = match angMom_b with - | (_,_,1) -> angMom_bz, 2 - | (_,1,_) -> angMom_by, 1 - | (_,_,_) -> angMom_bx, 0 + | (_,_,1) -> abz, 2 + | (_,1,_) -> aby, 1 + | (_,_,_) -> abx, 0 in if (bxyz < 1) then empty 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) + | 0 -> (aax+1,aay,aaz),(abx-1,aby,abz) + | 1 -> (aax,aay+1,aaz),(abx,aby-1,abz) + | _ -> (aax,aay,aaz+1),(abx,aby,abz-1) in let h1 = @@ -250,20 +252,20 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) else hrr0_v 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 (acx, acy, acz) = angMom_c + and (adx, ady, adz) = 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 + | (_,0,0) -> (acx+1, acy, acz), (adx-1, ady, adz), 0 + | (_,_,0) -> (acx, acy+1, acz), (adx, ady-1, adz), 1 + | _ -> (acx, acy, acz+1), (adx, ady, adz-1), 2 in let h1 = hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) and h2 = hrr_v angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) in - Array.mapi (fun k center_cd -> h1.(k) +. h2.(k) *. (Coordinate.coord center_cd xyz)) center_cd + Array.init ncoef (fun k -> h1.(k) +. h2.(k) *. (Coordinate.coord center_cd.(k) xyz)) in hrr_v (angMom_a.(0),angMom_a.(1),angMom_a.(2))