diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 9d6d6d7..bb6f559 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -12,12 +12,8 @@ open Util 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 + 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) *. @@ -34,8 +30,14 @@ type n_cls = { n : int ; cls : Z.t array } (** Write all integrals to a file with the convention *) let to_file ~filename basis = + 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 - let key_array = Array.make 12 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 @@ -46,48 +48,36 @@ let to_file ~filename basis = contracted_class basis.(i) basis.(j) basis.(k) basis.(l) in - for i_c = 0 to (Array.length basis.(i).Contracted_shell.powers) - 1 do - let x = Zkey.(to_int_array Kind_3 basis.(i).Contracted_shell.powers.(i_c)) in - key_array.(0) <- x.(0); - key_array.(1) <- x.(1); - key_array.(2) <- x.(2); - for j_c = 0 to (Array.length basis.(j).Contracted_shell.powers) - 1 do - let x = Zkey.(to_int_array Kind_3 basis.(j).Contracted_shell.powers.(j_c)) in - key_array.(3) <- x.(0); - key_array.(4) <- x.(1); - key_array.(5) <- x.(2); - for k_c = 0 to (Array.length basis.(k).Contracted_shell.powers) - 1 do - let x = Zkey.(to_int_array Kind_3 basis.(k).Contracted_shell.powers.(k_c)) in - key_array.(6) <- x.(0); - key_array.(7) <- x.(1); - key_array.(8) <- x.(2); - for l_c = 0 to (Array.length basis.(l).Contracted_shell.powers) - 1 do - let x = Zkey.(to_int_array Kind_3 basis.(l).Contracted_shell.powers.(l_c)) in - key_array.( 9) <- x.(0); - key_array.(10) <- x.(1); - key_array.(11) <- x.(2); - let key = - Zkey.(of_int_array Kind_12 key_array) - in - let value = - Zmap.find cls key - in - if (abs_float value > cutoff) then - Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" - (basis.(i).Contracted_shell.indice + i_c + 1) - (basis.(k).Contracted_shell.indice + k_c + 1) - (basis.(j).Contracted_shell.indice + j_c + 1) - (basis.(l).Contracted_shell.indice + l_c + 1) - value - done - done - done - done; + (* 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 = + Zkey.of_int_tuple (Zkey.Twelve (xi,xj,xk,xl)) + in + let value = + Zmap.find cls key + in + if (abs_float value > cutoff) 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; done; - done - ; + done; close_out oc (*