From 6d78d155f4a01623905ecc9086e32ce4ab00a0fb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 8 Mar 2019 17:47:58 +0100 Subject: [PATCH] Added functor for Two-electron integrals --- Basis/ERI.ml | 349 +++------------------------------ Basis/ERI.mli | 19 -- Basis/ERI_parallel.ml | 267 ------------------------- Basis/TwoElectronIntegrals.ml | 314 +++++++++++++++++++++++++++++ Basis/TwoElectronIntegrals.mli | 47 +++++ Makefile.include | 2 +- 6 files changed, 386 insertions(+), 612 deletions(-) delete mode 100644 Basis/ERI.mli delete mode 100644 Basis/ERI_parallel.ml create mode 100644 Basis/TwoElectronIntegrals.ml create mode 100644 Basis/TwoElectronIntegrals.mli diff --git a/Basis/ERI.ml b/Basis/ERI.ml index a9b3aae..8e7638b 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -1,336 +1,35 @@ (** Electron-electron repulsion integrals *) -open Util open Constants +open Util -include FourIdxStorage +module Zm = struct + let name = "Electron repulsion integrals" -module Am = AngularMomentum -module As = AtomicShell -module Asp = AtomicShellPair -module Bs = Basis -module Cs = ContractedShell -module Csp = ContractedShellPair -module Cspc = ContractedShellPairCouple -module Fis = FourIdxStorage - - -let cutoff = integrals_cutoff - - -(** (00|00)^m : Fundamental electron repulsion integral - $ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $ - - 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 = - assert (expo_pq_inv <> 0.); - let norm_pq_sq = - if norm_pq_sq > integrals_cutoff then norm_pq_sq else 0. - in - let exp_pq = 1. /. expo_pq_inv in - let t = norm_pq_sq *. exp_pq in - let f = two_over_sq_pi *. (sqrt exp_pq) in - let result = boys_function ~maxm t in - let rec aux accu k = function - | 0 -> result.(k) <- result.(k) *. accu - | l -> - begin - result.(k) <- result.(k) *. accu; - let new_accu = -. accu *. exp_pq in - aux new_accu (k+1) (l-1) - end - in - aux f 0 maxm; - result - - - - - -let class_of_contracted_shell_pair_couple shell_pair_couple = - let shell_p = Cspc.shell_pair_p shell_pair_couple - and shell_q = Cspc.shell_pair_q shell_pair_couple - in - if Array.length (Csp.shell_pairs shell_p) + (Array.length (Csp.shell_pairs shell_q)) < 4 then - TwoElectronRR.contracted_class_shell_pair_couple ~zero_m shell_pair_couple - else - TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m shell_p shell_q - - - - - -let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs = - List.map (fun pair -> - match Cspc.make ~cutoff pair pair with - | Some cspc -> - let cls = class_of_contracted_shell_pair_couple cspc in - (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) - (* TODO \sum_k |coef_k * integral_k| *) - | None -> (pair, -1.) - ) shell_pairs - |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) - |> List.map fst - - -(* TODO -let filter_contracted_shell_pair_couples ?(cutoff=integrals_cutoff) shell_pair_couples = - List.map (fun pair -> - let cls = - class_of_contracted_shell_pairs pair pair - in - (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) - ) shell_pairs - |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) - |> List.map fst - -*) - - -let store_class ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls = - let to_powers x = - let open Zkey in - match to_powers x with - | Three x -> x - | _ -> assert false - in - - let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple - and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple - in - - Array.iteri (fun i_c powers_i -> - let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in - let xi = to_powers powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in - let xj = to_powers powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in - let xk = to_powers powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in - let xl = to_powers powers_l in - let key = Zkey.of_powers_twelve xi xj xk xl in - let value = Zmap.find cls key in - set_chem data i_c j_c k_c l_c value - ) (Cs.zkey_array (Csp.shell_b shell_q)) - ) (Cs.zkey_array (Csp.shell_a shell_q)) - ) (Cs.zkey_array (Csp.shell_b shell_p)) - ) (Cs.zkey_array (Csp.shell_a shell_p)) - - - - - - -let of_basis_serial basis = - - let n = Bs.size basis - and shell = Bs.contracted_shells basis - in - - let eri_array = - Fis.create ~size:n `Dense -(* - Fis.create ~size:n `Sparse -*) - in - - let t0 = Unix.gettimeofday () in - - let shell_pairs = - Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs ~cutoff - in - - Printf.printf "%d significant shell pairs computed in %f seconds\n" - (List.length shell_pairs) (Unix.gettimeofday () -. t0); - - - let t0 = Unix.gettimeofday () in - let ishell = ref 0 in - - List.iter (fun shell_p -> - let () = - if (Cs.index (Csp.shell_a shell_p) > !ishell) then - (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) - in - - let sp = - Csp.shell_pairs shell_p - in - - try - List.iter (fun shell_q -> - let () = - if Cs.index (Csp.shell_a shell_q) > - Cs.index (Csp.shell_a shell_p) then - raise Exit - in - let sq = Csp.shell_pairs shell_q in - let cspc = - if Array.length sp < Array.length sq then - Cspc.make ~cutoff shell_p shell_q - else - Cspc.make ~cutoff shell_q shell_p - in - - match cspc with - | Some cspc -> - let cls = - class_of_contracted_shell_pair_couple cspc - in - store_class ~cutoff eri_array cspc cls - | None -> () - ) shell_pairs - with Exit -> () - ) shell_pairs ; - Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0); - eri_array - - - - - - - -(* Parallel functions *) - - - -let of_basis_parallel basis = - - let n = Bs.size basis - and shell = Bs.contracted_shells basis - in - - let store_class_parallel - ?(cutoff=integrals_cutoff) contracted_shell_pair_couple cls = - let to_powers x = - let open Zkey in - match to_powers x with - | Three x -> x - | _ -> assert false + let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = + assert (expo_pq_inv <> 0.); + let norm_pq_sq = + if norm_pq_sq > integrals_cutoff then norm_pq_sq else 0. in - - let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple - and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple + let exp_pq = 1. /. expo_pq_inv in + let t = norm_pq_sq *. exp_pq in + let f = two_over_sq_pi *. (sqrt exp_pq) in + let result = boys_function ~maxm t in + let rec aux accu k = function + | 0 -> result.(k) <- result.(k) *. accu + | l -> + begin + result.(k) <- result.(k) *. accu; + let new_accu = -. accu *. exp_pq in + aux new_accu (k+1) (l-1) + end in + aux f 0 maxm; + result - let result = ref [] in - Array.iteri (fun i_c powers_i -> - let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in - let xi = to_powers powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in - let xj = to_powers powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in - let xk = to_powers powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in - let xl = to_powers powers_l in - let key = Zkey.of_powers_twelve xi xj xk xl in - let value = Zmap.find cls key in - result := (i_c, j_c, k_c, l_c, value) :: !result - ) (Cs.zkey_array (Csp.shell_b shell_q)) - ) (Cs.zkey_array (Csp.shell_a shell_q)) - ) (Cs.zkey_array (Csp.shell_b shell_p)) - ) (Cs.zkey_array (Csp.shell_a shell_p)); - !result - in - - - - let t0 = Unix.gettimeofday () in - - let shell_pairs = - Csp.of_contracted_shell_array shell - |> filter_contracted_shell_pairs ~cutoff - in - - if Parallel.master then - Printf.printf "%d significant shell pairs computed in %f seconds\n" - (List.length shell_pairs) (Unix.gettimeofday () -. t0); - - - let t0 = Unix.gettimeofday () in - let ishell = ref max_int in - - let input_stream = Stream.of_list (List.rev shell_pairs) in - - let f shell_p = - let () = - if Parallel.rank < 2 && Cs.index (Csp.shell_a shell_p) < !ishell then - (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) - in - - let sp = - Csp.shell_pairs shell_p - in - - let result = ref [] in - try - List.iter (fun shell_q -> - let () = - if Cs.index (Csp.shell_a shell_q) > - Cs.index (Csp.shell_a shell_p) then - raise Exit - in - let sq = Csp.shell_pairs shell_q in - let cspc = - if Array.length sp < Array.length sq then - Cspc.make ~cutoff shell_p shell_q - else - Cspc.make ~cutoff shell_q shell_p - in - - match cspc with - | Some cspc -> - let cls = - class_of_contracted_shell_pair_couple cspc - in - result := (store_class_parallel ~cutoff cspc cls) :: !result; - | None -> () - ) shell_pairs; - raise Exit - with Exit -> List.concat !result |> Array.of_list - in - - let eri_array = - if Parallel.master then - Fis.create ~size:n `Dense - else - Fis.create ~size:0 `Dense - in - Farm.run ~ordered:false ~f input_stream - |> Stream.iter (fun l -> - Array.iter (fun (i_c,j_c,k_c,l_c,value) -> - set_chem eri_array i_c j_c k_c l_c value) l); - - if Parallel.master then - Printf.printf "Computed ERIs in parallel in %f seconds\n%!" (Unix.gettimeofday () -. t0); - Parallel.broadcast (lazy eri_array) - - - -let of_basis = - match Parallel.size with - | 1 -> of_basis_serial - | _ -> of_basis_parallel - - - - - - +end +module M = TwoElectronIntegrals.Make(Zm) +include M diff --git a/Basis/ERI.mli b/Basis/ERI.mli deleted file mode 100644 index 67206cf..0000000 --- a/Basis/ERI.mli +++ /dev/null @@ -1,19 +0,0 @@ -(** Electron repulsion intergals. *) - -include module type of FourIdxStorage - -(* -val filter_atomic_shell_pairs : ?cutoff:float -> AtomicShellPair.t list -> AtomicShellPair.t list -(** Uses Schwartz screening on atomic shell pairs. *) -*) - -val filter_contracted_shell_pairs : ?cutoff:float -> ContractedShellPair.t list -> ContractedShellPair.t list -(** Uses Schwartz screening on contracted shell pairs. *) - -val class_of_contracted_shell_pair_couple : ContractedShellPairCouple.t -> float Zmap.t -(** Computes all the ERI of the class built from a couple of contracted shell pairs. *) - -val of_basis : Basis.t -> t -(** Compute all ERI's for a given {!Basis.t}. *) - - diff --git a/Basis/ERI_parallel.ml b/Basis/ERI_parallel.ml deleted file mode 100644 index 5418507..0000000 --- a/Basis/ERI_parallel.ml +++ /dev/null @@ -1,267 +0,0 @@ -(** Electron-electron repulsion integrals *) - -open Util -open Constants -open Bigarray - -type t = (float, float32_elt, fortran_layout) Genarray.t - -(* Input type *) -type input_data = - { - i : int; - j : int; - shell_pairs : ContractedShellPair.t array array ; - schwartz : (float Zmap.t * float) array array; - cutoff : float - } - -(* Output type *) - -type output_integral = (* *) - { - i1 : int ; (* Function i for electron 1 *) - j2 : int ; (* Function j for electron 2 *) - k1 : int ; (* Function k for electron 1 *) - l2 : int ; (* Function l for electron 2 *) - swap : bool ; (* If true, Compute (kl|ij) instead of (ij|kl) *) - cls : float Zmap.t; - } - -type output_data = output_integral list - - -(** (00|00)^m : Fundamental electron repulsion integral - $ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $ - - 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 - let f = two_over_sq_pi *. (sqrt exp_pq) in - let result = boys_function ~maxm t in - let rec aux accu k = function - | 0 -> result.(k) <- result.(k) *. accu - | l -> - begin - result.(k) <- result.(k) *. accu; - let new_accu = -. accu *. exp_pq in - aux new_accu (k+1) (l-1) - end - in - aux f 0 maxm; - result - - -(** Compute all the integrals of a contracted class *) -let contracted_class_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - TwoElectronRR.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q - -(** Compute all the integrals of a contracted class *) -let contracted_class_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q - - - -(** Compute all the integrals of an atomic shell pair couple *) -let contracted_class_atomic_shell_pairs ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - TwoElectronRR.contracted_class_atomic_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q - -(** Compute all the integrals of a contracted class *) -let contracted_class_atomic_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - TwoElectronRR.contracted_class_atomic_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q - - - -(** Creates the input data structure from the basis set. *) -let make_input basis = - - let shell = Basis.contracted_shells basis in - - (* Pre-compute all shell pairs *) - let shell_pairs = - ContractedShellPair.shell_pairs shell - in - - (* Pre-compute diagonal integrals for Schwartz screening *) - let t0 = Unix.gettimeofday () in - - let schwartz = - Array.map (fun pair_array -> - Array.map (fun pair -> - let cls = contracted_class_shell_pairs pair pair in - let m = Zmap.fold (fun _ value accu -> - max (abs_float value) accu) cls 0. - in (cls, m) - ) pair_array - ) shell_pairs - in - - (* Count number of significant shell pairs. *) - let icount = ref 0 in - for i=0 to (Array.length shell) - 1 do - print_int (Contracted_shell.index shell.(i)) ; 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; - done; - done; - Printf.printf "%d shell pairs computed in %f seconds\n" !icount (Unix.gettimeofday () -. t0); - - List.init (Array.length shell) (fun i -> - List.init (i+1) (fun j -> - { i ; j ; shell_pairs ; schwartz ; cutoff } ) ) - |> List.fold_left List.rev_append [] - - - -exception NullIntegral - -let slave_job { i ; j ; shell_pairs ; schwartz ; cutoff} = - - let shell_p = shell_pairs.(i).(j) in - let schwartz_p, schwartz_p_max = schwartz.(i).(j) in - if schwartz_p_max < cutoff then [] else - let schwartz_cutoff = cutoff *. cutoff in - let sp = shell_p.ContractedShellPair.shell_pairs in - - let f k l = - let schwartz_q, schwartz_q_max = schwartz.(k).(l) in - if schwartz_p_max *. schwartz_q_max < schwartz_cutoff then - raise NullIntegral; - - let shell_q = shell_pairs.(k).(l) in - let sq = shell_q.ContractedShellPair.shell_pairs in - - let swap = Array.length sp > Array.length sq in - - (* Compute all the integrals of the class *) - let cls = - if swap then - if Array.length sp + Array.length sq = 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 - else - if Array.length sp + Array.length sq = 2 then - contracted_class_shell_pairs ~schwartz_p ~schwartz_q shell_p shell_q - else - contracted_class_shell_pairs_vec ~schwartz_p ~schwartz_q shell_p shell_q - in - { i1=i ; j2=k ; k1=j ; l2=l ; swap ; cls } - in - let rec loop accu k l = - match k, l with - | -1, -1 -> accu - | k, -1 -> loop accu (k-1) (k-1) - | k, l -> - let new_accu = - let _, schwartz_q_max = schwartz.(k).(l) in - if schwartz_p_max *. schwartz_q_max > schwartz_cutoff then - f k l :: accu - else accu - in - loop new_accu k (l-1) - in - loop [] i i - - -let of_basis basis = - - let shell = Basis.contracted_shells basis in - - (* 4D data initialization *) - let eri_array = - let n = Basis.size basis in - Genarray.create Float32 fortran_layout [| n ; n ; n ; n|] - in - Genarray.fill eri_array 0.; - - (* Compute ERIs *) - - let to_int_tuple x = - let open Zkey in - match to_int_tuple Kind_3 x with - | Three x -> x - | _ -> assert false - in - - let t0 = Unix.gettimeofday () in - let inn = ref 0 in - let out = ref 0 in - - make_input basis - |> List.map slave_job - |> List.iter (fun output_ij -> - List.iter (fun { i1=i ; j2=k ; k1=j ; l2=l ; swap ; cls } -> - - Array.iteri (fun i_c powers_i -> - let i_c = (Contracted_shell.index shell.(i)) + i_c + 1 in - let xi = to_int_tuple powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = (Contracted_shell.index shell.(j)) + j_c + 1 in - let xj = to_int_tuple powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = (Contracted_shell.index shell.(k)) + k_c + 1 in - let xk = to_int_tuple powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = (Contracted_shell.index shell.(l)) + 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 - eri_array.{i_c,k_c,j_c,l_c} <- value; - eri_array.{j_c,k_c,i_c,l_c} <- value; - eri_array.{i_c,l_c,j_c,k_c} <- value; - eri_array.{j_c,l_c,i_c,k_c} <- value; - eri_array.{k_c,i_c,l_c,j_c} <- value; - eri_array.{k_c,j_c,l_c,i_c} <- value; - eri_array.{l_c,i_c,k_c,j_c} <- value; - eri_array.{l_c,j_c,k_c,i_c} <- value; - if (abs_float value > cutoff) then - (inn := !inn + 1; - ) - else - out := !out + 1; - ) (Contracted_shell.powers shell.(l)) - ) (Contracted_shell.powers shell.(k)) - ) (Contracted_shell.powers shell.(j)) - ) (Contracted_shell.powers shell.(i)); - ) output_ij - ); - Printf.printf "In: %d Out:%d\n" !inn !out ; - Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0); - eri_array - - -(** Write all integrals to a file with the convention *) -let to_file ~filename eri_array = - let oc = open_out filename in - (* Print ERIs *) - for l_c=1 to (Genarray.nth_dim eri_array 3) do - for k_c=1 to l_c do - for j_c=1 to l_c do - for i_c=1 to k_c do - let value = eri_array.{i_c,j_c,k_c,l_c} in - if (abs_float value > cutoff) then - Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i_c j_c k_c l_c value; - done; - done; - done; - done; - close_out oc - - diff --git a/Basis/TwoElectronIntegrals.ml b/Basis/TwoElectronIntegrals.ml new file mode 100644 index 0000000..34dcf7b --- /dev/null +++ b/Basis/TwoElectronIntegrals.ml @@ -0,0 +1,314 @@ +(** Two electron integral functor, parameterized by the zero_m function. +*) + +open Constants +let cutoff = integrals_cutoff + +module Am = AngularMomentum +module As = AtomicShell +module Asp = AtomicShellPair +module Bs = Basis +module Cs = ContractedShell +module Csp = ContractedShellPair +module Cspc = ContractedShellPairCouple +module Fis = FourIdxStorage + + +module type Zero_mType = + sig + val name : string + val zero_m : maxm:int -> expo_pq_inv:float -> norm_pq_sq:float -> float array + end + + +module Make(Zero_m : Zero_mType) = struct + + include FourIdxStorage + + let zero_m = Zero_m.zero_m + + let class_of_contracted_shell_pair_couple shell_pair_couple = + let shell_p = Cspc.shell_pair_p shell_pair_couple + and shell_q = Cspc.shell_pair_q shell_pair_couple + in + if Array.length (Csp.shell_pairs shell_p) + + (Array.length (Csp.shell_pairs shell_q)) < 4 then + TwoElectronRR.contracted_class_shell_pair_couple + ~zero_m shell_pair_couple + else + TwoElectronRRVectorized.contracted_class_shell_pairs + ~zero_m shell_p shell_q + + + + let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs = + List.map (fun pair -> + match Cspc.make ~cutoff pair pair with + | Some cspc -> + let cls = class_of_contracted_shell_pair_couple cspc in + (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) + (* TODO \sum_k |coef_k * integral_k| *) + | None -> (pair, -1.) + ) shell_pairs + |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) + |> List.map fst + + + (* TODO + let filter_contracted_shell_pair_couples + ?(cutoff=integrals_cutoff) shell_pair_couples = + List.map (fun pair -> + let cls = + class_of_contracted_shell_pairs pair pair + in + (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) + ) shell_pairs + |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) + |> List.map fst + + *) + + + let store_class ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls = + let to_powers x = + let open Zkey in + match to_powers x with + | Three x -> x + | _ -> assert false + in + + let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple + and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple + in + + Array.iteri (fun i_c powers_i -> + let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in + let xi = to_powers powers_i in + Array.iteri (fun j_c powers_j -> + let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in + let xj = to_powers powers_j in + Array.iteri (fun k_c powers_k -> + let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in + let xk = to_powers powers_k in + Array.iteri (fun l_c powers_l -> + let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in + let xl = to_powers powers_l in + let key = Zkey.of_powers_twelve xi xj xk xl in + let value = Zmap.find cls key in + set_chem data i_c j_c k_c l_c value + ) (Cs.zkey_array (Csp.shell_b shell_q)) + ) (Cs.zkey_array (Csp.shell_a shell_q)) + ) (Cs.zkey_array (Csp.shell_b shell_p)) + ) (Cs.zkey_array (Csp.shell_a shell_p)) + + + + + + + let of_basis_serial basis = + + let n = Bs.size basis + and shell = Bs.contracted_shells basis + in + + let eri_array = + Fis.create ~size:n `Dense +(* + Fis.create ~size:n `Sparse +*) + in + + let t0 = Unix.gettimeofday () in + + let shell_pairs = + Csp.of_contracted_shell_array shell + |> filter_contracted_shell_pairs ~cutoff + in + + Printf.printf "%d significant shell pairs computed in %f seconds\n" + (List.length shell_pairs) (Unix.gettimeofday () -. t0); + + + let t0 = Unix.gettimeofday () in + let ishell = ref 0 in + + List.iter (fun shell_p -> + let () = + if (Cs.index (Csp.shell_a shell_p) > !ishell) then + (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) + in + + let sp = + Csp.shell_pairs shell_p + in + + try + List.iter (fun shell_q -> + let () = + if Cs.index (Csp.shell_a shell_q) > + Cs.index (Csp.shell_a shell_p) then + raise Exit + in + let sq = Csp.shell_pairs shell_q in + let cspc = + if Array.length sp < Array.length sq then + Cspc.make ~cutoff shell_p shell_q + else + Cspc.make ~cutoff shell_q shell_p + in + + match cspc with + | Some cspc -> + let cls = + class_of_contracted_shell_pair_couple cspc + in + store_class ~cutoff eri_array cspc cls + | None -> () + ) shell_pairs + with Exit -> () + ) shell_pairs ; + Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0); + eri_array + + + + + + + + (* Parallel functions *) + + + + let of_basis_parallel basis = + + let n = Bs.size basis + and shell = Bs.contracted_shells basis + in + + let store_class_parallel + ?(cutoff=integrals_cutoff) contracted_shell_pair_couple cls = + let to_powers x = + let open Zkey in + match to_powers x with + | Three x -> x + | _ -> assert false + in + + let shell_p = Cspc.shell_pair_p contracted_shell_pair_couple + and shell_q = Cspc.shell_pair_q contracted_shell_pair_couple + in + + let result = ref [] in + Array.iteri (fun i_c powers_i -> + let i_c = Cs.index (Csp.shell_a shell_p) + i_c + 1 in + let xi = to_powers powers_i in + Array.iteri (fun j_c powers_j -> + let j_c = Cs.index (Csp.shell_b shell_p) + j_c + 1 in + let xj = to_powers powers_j in + Array.iteri (fun k_c powers_k -> + let k_c = Cs.index (Csp.shell_a shell_q) + k_c + 1 in + let xk = to_powers powers_k in + Array.iteri (fun l_c powers_l -> + let l_c = Cs.index (Csp.shell_b shell_q) + l_c + 1 in + let xl = to_powers powers_l in + let key = Zkey.of_powers_twelve xi xj xk xl in + let value = Zmap.find cls key in + result := (i_c, j_c, k_c, l_c, value) :: !result + ) (Cs.zkey_array (Csp.shell_b shell_q)) + ) (Cs.zkey_array (Csp.shell_a shell_q)) + ) (Cs.zkey_array (Csp.shell_b shell_p)) + ) (Cs.zkey_array (Csp.shell_a shell_p)); + !result + in + + + + let t0 = Unix.gettimeofday () in + + let shell_pairs = + Csp.of_contracted_shell_array shell + |> filter_contracted_shell_pairs ~cutoff + in + + if Parallel.master then + Printf.printf "%d significant shell pairs computed in %f seconds\n" + (List.length shell_pairs) (Unix.gettimeofday () -. t0); + + + let t0 = Unix.gettimeofday () in + let ishell = ref max_int in + + let input_stream = Stream.of_list (List.rev shell_pairs) in + + let f shell_p = + let () = + if Parallel.rank < 2 && Cs.index (Csp.shell_a shell_p) < !ishell then + (ishell := Cs.index (Csp.shell_a shell_p) ; print_int !ishell ; print_newline ()) + in + + let sp = + Csp.shell_pairs shell_p + in + + let result = ref [] in + try + List.iter (fun shell_q -> + let () = + if Cs.index (Csp.shell_a shell_q) > + Cs.index (Csp.shell_a shell_p) then + raise Exit + in + let sq = Csp.shell_pairs shell_q in + let cspc = + if Array.length sp < Array.length sq then + Cspc.make ~cutoff shell_p shell_q + else + Cspc.make ~cutoff shell_q shell_p + in + + match cspc with + | Some cspc -> + let cls = + class_of_contracted_shell_pair_couple cspc + in + result := (store_class_parallel ~cutoff cspc cls) :: !result; + | None -> () + ) shell_pairs; + raise Exit + with Exit -> List.concat !result |> Array.of_list + in + + let eri_array = + if Parallel.master then + Fis.create ~size:n `Dense + else + Fis.create ~size:0 `Dense + in + Farm.run ~ordered:false ~f input_stream + |> Stream.iter (fun l -> + Array.iter (fun (i_c,j_c,k_c,l_c,value) -> + set_chem eri_array i_c j_c k_c l_c value) l); + + if Parallel.master then + Printf.printf + "Computed %s Integrals in parallel in %f seconds\n%!" Zero_m.name (Unix.gettimeofday () -. t0); + Parallel.broadcast (lazy eri_array) + + + + let of_basis = + match Parallel.size with + | 1 -> of_basis_serial + | _ -> of_basis_parallel + + +end + + + + + + + diff --git a/Basis/TwoElectronIntegrals.mli b/Basis/TwoElectronIntegrals.mli new file mode 100644 index 0000000..db90cff --- /dev/null +++ b/Basis/TwoElectronIntegrals.mli @@ -0,0 +1,47 @@ +(** Two-electron integrals with an arbitrary operator, with a functorial interface + parameterized by the fundamental two-electron integrals. + + {% $(00|00)^m = \int \int \phi_p(r1) \hat{O} \phi_q(r2) dr_1 dr_2 $ %} : Fundamental two-electron integral + +*) + + +module type Zero_mType = + sig +val name : string +(** Name of the kind of integrals, for printing purposes. *) + +val zero_m : maxm:int -> expo_pq_inv:float -> norm_pq_sq:float -> float array +(** The returned float array contains all the {% $(00|00)^m$ %} values, where + [m] is the index of the array. + + - [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$ %} + +*) + + end + + +module Make : functor (Zero_m : Zero_mType) -> + sig + include module type of FourIdxStorage + + val class_of_contracted_shell_pair_couple : + ContractedShellPairCouple.t -> float Zmap.t + (** Computes all the ERI of the class built from a couple of + contracted shell pairs. *) + + val filter_contracted_shell_pairs : + ?cutoff:float -> + ContractedShellPair.t list -> ContractedShellPair.t list + (** Uses Schwartz screening on contracted shell pairs. *) + + val of_basis : Basis.t -> t + (** Compute all ERI's for a given {!Basis.t}. *) + + end + diff --git a/Makefile.include b/Makefile.include index 7063e7c..235218a 100644 --- a/Makefile.include +++ b/Makefile.include @@ -49,7 +49,7 @@ doc: QCaml.odocl $(OCAMLBUILD) -ocamlc ocamlcp $*.byte -use-ocamlfind $(PKGS) clean: - rm QCaml.odocl && $(OCAMLBUILD) -clean + rm -f QCaml.odocl && $(OCAMLBUILD) -clean debug: run_integrals.native ./debug.sh