From 925cdf693f8e06356470e70b800445d774fce2bc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 3 Feb 2018 23:26:20 +0100 Subject: [PATCH] Working on Nuclear integrals --- Basis/ERI.ml | 179 ++++++++++++++----------- Basis/NucInt.ml | 129 ++++++++++++++++++ Basis/OneElectronRR.ml | 290 ++++++++++++++++++++++------------------- Basis/Overlap.ml | 9 +- Basis/TwoElectronRR.ml | 1 - run_integrals.ml | 8 +- 6 files changed, 400 insertions(+), 216 deletions(-) create mode 100644 Basis/NucInt.ml diff --git a/Basis/ERI.ml b/Basis/ERI.ml index b25b24b..cf2b5a8 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -25,6 +25,10 @@ let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = (** Compute all the integrals of a contracted class *) +(* +let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t = + TwoElectronRRVectorized.contracted_class ~zero_m shell_a shell_b shell_c shell_d + *) let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t = TwoElectronRR.contracted_class ~zero_m shell_a shell_b shell_c shell_d @@ -37,6 +41,9 @@ let contracted_class_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float let cutoff2 = cutoff *. cutoff +(* +type n_cls = { n : int ; cls : Z.t array } +*) exception NullIntegral (* @@ -66,7 +73,10 @@ let to_file ~filename basis = (* Pre-compute all shell pairs *) let shell_pairs = Array.mapi (fun i shell_a -> Array.map (fun shell_b -> + (* Shell_pair.create_array shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis + *) + Shell_pair.create_array shell_a shell_b) (basis) ) basis in Printf.printf "%d shells\n" (Array.length basis); @@ -87,8 +97,8 @@ let to_file ~filename basis = print_int basis.(i).Contracted_shell.indice ; print_newline (); for j=0 to i do let schwartz_p, schwartz_p_max = schwartz.(i).(j) in - if (schwartz_p_max >= cutoff) then - icount := !icount + 1; + if (schwartz_p_max >= cutoff) then + icount := !icount + 1; done; done; Printf.printf "%d shell pairs computed in %f seconds\n" !icount (Unix.gettimeofday () -. t0); @@ -111,83 +121,111 @@ let to_file ~filename basis = 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 j=0 to (Array.length basis) - 1 do + (* let schwartz_p, schwartz_p_max = schwartz.(i).(j) in + *) try + (* if (schwartz_p_max < cutoff) then raise NullIntegral; + *) let shell_p = shell_pairs.(i).(j) in +(* for k=0 to i do for l=0 to k do +*) + for k=0 to (Array.length basis) - 1 do + for l=0 to (Array.length basis) - 1 do + (* let schwartz_q, schwartz_q_max = schwartz.(k).(l) in + *) try - if schwartz_p_max *. schwartz_q_max < cutoff2 then - raise NullIntegral; - let - shell_q = shell_pairs.(k).(l) - in + (* + if schwartz_p_max *. schwartz_q_max < cutoff2 then + raise NullIntegral; + *) + let + shell_q = shell_pairs.(k).(l) + in - let swap = - Array.length shell_q < Array.length shell_p - in + let swap = + Array.length shell_q < Array.length shell_p + in - (* Compute all the integrals of the class *) - let cls = - if swap then - if Array.length shell_p < 2 then - contracted_class_shell_pairs - ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p - shell_q shell_p - else - contracted_class_shell_pairs_vec - ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p - shell_q shell_p + (* Compute all the integrals of the class *) + let cls = + if swap then + if Array.length shell_p < 2 then + (* + contracted_class_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p + *) + contracted_class_shell_pairs shell_q shell_p else + (* + contracted_class_shell_pairs_vec ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p + *) + contracted_class_shell_pairs_vec shell_q shell_p + else if Array.length shell_q < 2 then - contracted_class_shell_pairs - ~schwartz_p ~schwartz_q - shell_p shell_q + (* + contracted_class_shell_pairs ~schwartz_p ~schwartz_q shell_p shell_q + *) + contracted_class_shell_pairs shell_p shell_q else - contracted_class_shell_pairs_vec - ~schwartz_p ~schwartz_q - shell_p shell_q - in + (* + contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q + *) + contracted_class_shell_pairs_vec shell_p shell_q + in - (* Write the data in the output file *) - Array.iteri (fun i_c powers_i -> - let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in - let xi = to_int_tuple powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in - let xj = to_int_tuple powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = basis.(k).Contracted_shell.indice + k_c + 1 in - let xk = to_int_tuple powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = basis.(l).Contracted_shell.indice + l_c + 1 in - let xl = to_int_tuple powers_l in - let key = - if swap then - Zkey.of_int_tuple (Zkey.Twelve (xk,xl,xi,xj)) - else - Zkey.of_int_tuple (Zkey.Twelve (xi,xj,xk,xl)) - in - let value = - Zmap.find cls key - in - if (abs_float value > cutoff) then - (inn := !inn + 1; - eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} <- value; - ) + (* Write the data in the output file *) + Array.iteri (fun i_c powers_i -> + let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in + let xi = to_int_tuple powers_i in + Array.iteri (fun j_c powers_j -> + let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in + let xj = to_int_tuple powers_j in + Array.iteri (fun k_c powers_k -> + let k_c = basis.(k).Contracted_shell.indice + k_c + 1 in + let xk = to_int_tuple powers_k in + Array.iteri (fun l_c powers_l -> + let l_c = basis.(l).Contracted_shell.indice + l_c + 1 in + let xl = to_int_tuple powers_l in + let key = + if swap then + Zkey.of_int_tuple (Zkey.Twelve (xk,xl,xi,xj)) else - out := !out + 1; - ) basis.(l).Contracted_shell.powers - ) basis.(k).Contracted_shell.powers - ) basis.(j).Contracted_shell.powers - ) basis.(i).Contracted_shell.powers; - with NullIntegral -> () + Zkey.of_int_tuple (Zkey.Twelve (xi,xj,xk,xl)) + in + let value = + Zmap.find cls key + in + eri_array.{(i_c-1),(k_c-1),(j_c-1),(l_c-1)} <- value; + (* + eri_array.{(j_c-1),(k_c-1),(i_c-1),(l_c-1)} <- value; + eri_array.{(i_c-1),(l_c-1),(j_c-1),(k_c-1)} <- value; + eri_array.{(j_c-1),(l_c-1),(i_c-1),(k_c-1)} <- value; + eri_array.{(k_c-1),(i_c-1),(l_c-1),(j_c-1)} <- value; + eri_array.{(k_c-1),(j_c-1),(l_c-1),(i_c-1)} <- value; + eri_array.{(l_c-1),(i_c-1),(k_c-1),(j_c-1)} <- value; + eri_array.{(l_c-1),(j_c-1),(k_c-1),(i_c-1)} <- value; + *) + if (abs_float value > cutoff) then + (inn := !inn + 1; + ) + else + out := !out + 1; + ) basis.(l).Contracted_shell.powers + ) basis.(k).Contracted_shell.powers + ) basis.(j).Contracted_shell.powers + ) basis.(i).Contracted_shell.powers; + with NullIntegral -> () done; done; with NullIntegral -> () @@ -198,15 +236,15 @@ let to_file ~filename basis = (* Print ERIs *) 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; + 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 (abs_float value > cutoff) 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 @@ -398,8 +436,3 @@ let xto_file ~filename basis = *) - - - - - diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml new file mode 100644 index 0000000..e551c8d --- /dev/null +++ b/Basis/NucInt.ml @@ -0,0 +1,129 @@ +(** Electron-nucleus repulsion integrals *) + +open Util +open Constants +open Bigarray + +(** (0|0)^m : Fundamental electron-nucleus repulsion integral + $ \int \phi_p(r1) 1/r_{C} dr_1 $ + + maxm : Maximum total angular momentum + expo_pq_inv : $1./p + 1./q$ where $p$ and $q$ are the exponents of + $\phi_p$ and $\phi_q$ + norm_pq_sq : square of the distance between the centers of $\phi_p$ + and $\phi_q$ +*) + +let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = + let exp_pq = 1. /. expo_pq_inv in + let t = norm_pq_sq *. exp_pq in + boys_function ~maxm t + |> Array.mapi (fun m fm -> + two_over_sq_pi *. (if m mod 2 = 0 then fm else -.fm) *. + (pow exp_pq m) *. (sqrt exp_pq) + ) + + +(** Compute all the integrals of a contracted class *) +let contracted_class_shell_pair shell_p nucl_coord : float Zmap.t = + OneElectronRR.contracted_class_shell_pair ~zero_m shell_p nucl_coord + + +let cutoff2 = cutoff *. cutoff +exception NullIntegral + +(* +(** Unique index for integral *) +let index i j k l = + let f i k = + let (p,r) = + if i <= k then (i,k) else (k,i) + in p+ (r*r-r)/2 + in + let p = f i k and q = f j l in + f p q +*) + + +(** Write all integrals to a file with the convention *) +let to_file ~filename basis nucl_coord = + let to_int_tuple x = + let open Zkey in + match to_int_tuple Kind_3 x with + | Three x -> x + | _ -> assert false + in + + let oc = open_out filename in + + (* Pre-compute all shell pairs *) + let shell_pairs = + Array.mapi (fun i shell_a -> Array.map (fun shell_b -> + Shell_pair.create_array shell_a shell_b) (Array.sub basis 0 (i+1)) ) basis + in + Printf.printf "%d shells\n" (Array.length basis); + + let eni_array = + let n = ref 0 in + for i=0 to (Array.length basis) - 1 do + n := !n + (Array.length (basis.(i).Contracted_shell.powers)) + done; + let n = !n in + Array2.create Float64 c_layout n n + in + Array2.fill eni_array 0.; + + (* Compute Integrals *) + let t0 = Unix.gettimeofday () in + + let inn = ref 0 and out = ref 0 in + + for i=0 to (Array.length basis) - 1 do + print_int basis.(i).Contracted_shell.indice ; print_newline (); + for j=0 to i do + let + shell_p = shell_pairs.(i).(j) + in + + (* Compute all the integrals of the class *) + let cls = + contracted_class_shell_pair shell_p nucl_coord + in + + (* Write the data in the output file *) + Array.iteri (fun i_c powers_i -> + let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in + let xi = to_int_tuple powers_i in + Array.iteri (fun j_c powers_j -> + let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in + let xj = to_int_tuple powers_j in + let key = + Zkey.of_int_tuple (Zkey.Six (xi,xj)) + in + let value = + Zmap.find cls key + in + if (abs_float value > cutoff) then + (inn := !inn + 1; + eni_array.{(i_c-1),(j_c-1)} <- value; + ) + else + out := !out + 1; + ) basis.(j).Contracted_shell.powers + ) basis.(i).Contracted_shell.powers; + done; + done; + + Printf.printf "Computed %d non-zero ENIs in %f seconds\n" !inn (Unix.gettimeofday () -. t0); + + (* Print ENIs *) + for i_c=1 to (Array2.dim1 eni_array) do + for j_c=1 to (Array2.dim2 eni_array) do + let value = eni_array.{(i_c-1),(j_c-1)} in + if (value <> 0.) then + Printf.fprintf oc " %5d %5d%20.15f\n" i_c j_c value; + done; + done; + Printf.printf "In: %d Out:%d\n" !inn !out ; + close_out oc + diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index 31aea67..26ae706 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -1,4 +1,7 @@ open Util +open Constants + +exception NullPair (** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *) let chop f g = @@ -9,7 +12,7 @@ let chop f g = (** Horizontal and Vertical Recurrence Relations (HVRR) *) let hvrr_one_e - m (angMom_a, angMom_b) (totAngMom_a, totAngMom_b) + (angMom_a, angMom_b) (totAngMom_a, totAngMom_b) (maxm, zero_m_array) (expo_inv_p) (center_ab, center_pa, center_pc) map = @@ -17,88 +20,99 @@ let hvrr_one_e 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 maxsze = maxm+1 in + let empty = Array.make maxsze 0. in (** Vertical recurrence relations *) - let rec vrr m angMom_a totAngMom_a = - - if angMom_a.(0) < 0 || angMom_a.(1) < 0 || angMom_a.(2) < 0 then 0. + let rec vrr angMom_a totAngMom_a = + let ax,ay,az = angMom_a in + if (ax < 0) || (ay < 0) || (az < 0) then + empty else match totAngMom_a with - | 0 -> zero_m_array.(m) + | 0 -> zero_m_array | _ -> + let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in - let key = [| angMom_a.(0)+1; angMom_a.(1)+1; angMom_a.(2)+1; |] - |> Zkey.(of_int_array ~kind:Kind_3) - in - - let (found, result) = - try (true, Zmap.find map.(m) key) with - | Not_found -> (false, - 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 - in - am.(xyz) <- am.(xyz) - 1; - amm.(xyz) <- amm.(xyz) - 2; - chop (Coordinate.coord center_pa xyz) - (fun () -> vrr m am (totAngMom_a-1)) - +. chop (0.5 *. (float_of_int am.(xyz)) *. expo_inv_p) - (fun () -> vrr m amm (totAngMom_a-2)) - -. chop ((Coordinate.coord center_pc xyz) *. expo_inv_p) - (fun () -> vrr (m+1) am (totAngMom_a-1)) - -. chop (0.5 *. (float_of_int am.(xyz)) *. expo_inv_p *. expo_inv_p) - (fun () -> vrr (m+1) amm (totAngMom_a-2)) - in - if not found then - Zmap.add map.(m) key result; - result + 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, 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 = Coordinate.coord center_pa xyz + and f2 = expo_inv_p *. (Coordinate.coord center_pc xyz) + in + if amxyz < 1 then + let v1 = + vrr am (totAngMom_a-1) + in + Array.init maxsze (fun m -> + if m = maxm then 0. else (f1 *. v1.(m) ) -. f2 *. v1.(m+1) ) + else + let v3 = + vrr amm (totAngMom_a-2) + in + let v1 = + vrr 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 key result; + result (** Horizontal recurrence relations *) and hrr angMom_a angMom_b totAngMom_a totAngMom_b = - if angMom_b.(0) < 0 || angMom_b.(1) < 0 || angMom_b.(2) < 0 then 0. + let bx,by,bz = angMom_b in + if (bx < 0) || (by < 0) || (bz < 0) then 0. else match totAngMom_b with - | 0 -> vrr 0 angMom_a + | 0 -> (vrr angMom_a totAngMom_a).(0) | _ -> - - let key = [| angMom_a.(0)+1; angMom_a.(1)+1; angMom_a.(2)+1; - angMom_b.(0)+1; angMom_b.(1)+1; angMom_b.(2)+1 |] - |> Zkey.(of_int_array ~kind:Kind_6) + 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 - - let (found, result) = - try (true, Zmap.find map.(m) key) with - | Not_found -> (false, - begin - 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 - in - ap.(xyz) <- ap.(xyz) + 1; - bm.(xyz) <- bm.(xyz) - 1; - hrr ap bm (totAngMom_a+1) (totAngMom_b-1) - +. chop (Coordinate.coord center_ab xyz) (fun () -> - hrr angMom_a bm totAngMom_a (totAngMom_b-1) ) - end) - in - if not found then - Zmap.add map.(m) key result; - result + 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 = + hrr ap bm (totAngMom_a+1) (totAngMom_b-1) + in + let f2 = + (Coordinate.coord center_ab xyz) + in + if (abs_float f2 < cutoff) then h1 else + let h2 = + hrr angMom_a bm totAngMom_a (totAngMom_b-1) + in + h1 +. f2 *. h2 in - hrr m 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)) + totAngMom_a totAngMom_b @@ -109,10 +123,12 @@ let hvrr_one_e (** Computes all the one-electron integrals of the contracted shell pair *) -let contracted_class_nuc ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = +let contracted_class_shell_pair ~zero_m shell_p nucl_coord : float Zmap.t = - let shell_p = Shell_pair.create_array shell_a shell_b - and maxm = + let shell_a = shell_p.(0).Shell_pair.shell_a + and shell_b = shell_p.(0).Shell_pair.shell_b + in + let maxm = let open Angular_momentum in (to_int @@ Contracted_shell.totAngMom shell_a) + (to_int @@ Contracted_shell.totAngMom shell_b) in @@ -127,80 +143,80 @@ let contracted_class_nuc ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t let contracted_class = Array.make (Array.length class_indices) 0.; in - () (* Compute all integrals in the shell for each pair of significant shell pairs *) for ab=0 to (Array.length shell_p - 1) do - let coef_prod = - shell_p.(ab).Shell_pair.coef - in - (** Screening on thr product of coefficients *) - if (abs_float coef_prod) > 1.e-4*.cutoff then - begin + try + begin + let coef_prod = shell_p.(ab).Shell_pair.coef in + let norm_coef_scale_p = shell_p.(ab).Shell_pair.norm_coef_scale in - let expo_pq_inv = - shell_p.(ab).Shell_pair.expo_inv + (** Screening on the product of coefficients *) + if (abs_float coef_prod) < 1.e-4*.cutoff then + raise NullPair; + + + let expo_pq_inv = + shell_p.(ab).Shell_pair.expo_inv + in + + let center_ab = + shell_p.(ab).Shell_pair.center_ab + in + let center_pa = + Coordinate.(center_ab |- shell_a.Contracted_shell.center) + in + + for c=0 to Array.length nucl_coord - 1 do + let center_pc = + Coordinate.(shell_p.(ab).Shell_pair.center |- nucl_coord.(c) ) + in + let norm_pq_sq = + Coordinate.dot center_pc center_pc + in + + let zero_m_array = + zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + in + match Contracted_shell.(totAngMom shell_a, totAngMom shell_b) with + | Angular_momentum.(S,S) -> + let integral = + zero_m_array.(0) in - let center_ab = - Coordinate.shell_p.(ab).Shell_pair.center_ab - in - let norm_pq_sq = - Coordinate.dot center_pq center_pq - in - - let zero_m_array = - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq - in - - match Contracted_shell.(totAngMom shell_a, totAngMom shell_b, - totAngMom shell_c, totAngMom shell_d) with - | Angular_momentum.(S,S,S,S) -> Array.iteri (fun i key -> - let coef_prod = - shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef - in - let integral = - zero_m_array.(0) - in - contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral - ) class_indices - | _ -> - let d = shell_q.(cd).Shell_pair.j in - let map = Array.init maxm (fun _ -> Zmap.create (Array.length class_indices)) in - (* Compute the integral class from the primitive shell quartet *) - Array.iteri (fun i key -> - let (angMomA,angMomB,angMomC,angMomD) = - let a = Zkey.to_int_array Zkey.Kind_12 key in - ( [| a.(0) ; a.(1) ; a.(2) |], - [| a.(3) ; a.(4) ; a.(5) |], - [| a.(6) ; a.(7) ; a.(8) |], - [| a.(9) ; a.(10) ; a.(11) |] ) - in - let norm = - shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD - in - let integral = chop norm (fun () -> - ghvrr 0 (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) - (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d) - (shell_p.(ab).Shell_pair.expo_inv, shell_q.(cd).Shell_pair.expo_inv) - (shell_p.(ab).Shell_pair.center_ab, shell_q.(cd).Shell_pair.center_ab, center_pq) - map ) - in - contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral - ) class_indices - end - done - done; - let result = - Zmap.create (Array.length contracted_class) - in - Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; - result - -*) + contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral + | _ -> + let map = Zmap.create (2*maxm) in + let norm_coef_scale = norm_coef_scale_p in + (* Compute the integral class from the primitive shell quartet *) + class_indices + |> Array.iteri (fun i key -> + let (angMomA,angMomB) = + let a = Zkey.to_int_array Zkey.Kind_6 key in + ( [| a.(0) ; a.(1) ; a.(2) |], + [| a.(3) ; a.(4) ; a.(5) |]) + in + let norm = norm_coef_scale.(i) in + 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) shell_p.(ab).Shell_pair.expo_inv + (center_ab, center_pa, center_pc) + map + in + contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral + ) + done + end + with NullPair -> () +done; +let result = + Zmap.create (Array.length contracted_class) +in +Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; +result + diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index 38dc425..276aa3c 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -1,4 +1,5 @@ open Util +open Constants (** Computes all the overlap integrals of the contracted shell pair *) let contracted_class shell_a shell_b : float Zmap.t = @@ -36,6 +37,9 @@ let contracted_class shell_a shell_b : float Zmap.t = let expo_inv = shell_p.(ab).Shell_pair.expo_inv in + let norm_coef_scale = + shell_p.(ab).Shell_pair.norm_coef_scale + in Array.iteri (fun i key -> let (angMomA,angMomB) = @@ -43,15 +47,14 @@ let contracted_class shell_a shell_b : float Zmap.t = ( [| a.(0) ; a.(1) ; a.(2) |], [| a.(3) ; a.(4) ; a.(5) |] ) in - let norm = - shell_p.(ab).Shell_pair.norm_fun angMomA angMomB - in let f k = Overlap_primitives.hvrr (angMomA.(k), angMomB.(k)) expo_inv (Coordinate.coord center_ab k, Coordinate.coord center_a k) in + let norm = norm_coef_scale.(i) in + let coef_prod = coef_prod *. norm in let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral ) class_indices diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index a5bc679..bdadb85 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -50,7 +50,6 @@ let hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) match totAngMom_a with | 0 -> zero_m_array | _ -> - let maxsze = maxm+1 in let key = Zkey.of_int_tuple (Zkey.Three angMom_a) in try Zmap.find map_1d key with diff --git a/run_integrals.ml b/run_integrals.ml index 96c1950..75641cc 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -27,10 +27,14 @@ let run ~out = let basis = Lazy.force Basis.basis in print_endline @@ Basis.to_string basis; - ERI.to_file ~filename:(out_file^".eri") basis -(* + (* Overlap.to_file ~filename:(out_file^".overlap") basis + ERI.to_file ~filename:(out_file^".eri") basis *) + let nucl_coord = + Array.map (fun (_,c) -> c) nuclei + in + NucInt.to_file ~filename:(out_file^".nuc") basis nucl_coord let () =