From 0face5817ec5079edd9c79c5eb23535e0b3822f6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 22 Mar 2018 00:29:14 +0100 Subject: [PATCH] More functional --- Basis/AtomicShellPair.ml | 8 +- Basis/AtomicShellPairCouple.ml | 5 +- Basis/ContractedShell.ml | 2 + Basis/ContractedShell.mli | 3 + Basis/ContractedShellPair.ml | 16 ++- Basis/ContractedShellPair.mli | 16 +-- Basis/ContractedShellPairCouple.ml | 5 +- Basis/ERI.ml | 205 ++++++++++++++++------------- Basis/PrimitiveShell.ml | 2 + Basis/PrimitiveShell.mli | 3 + Basis/PrimitiveShellPair.ml | 78 +++++------ Basis/PrimitiveShellPair.mli | 3 +- Utils/Util.ml | 9 ++ Utils/Util.mli | 4 + 14 files changed, 203 insertions(+), 156 deletions(-) diff --git a/Basis/AtomicShellPair.ml b/Basis/AtomicShellPair.ml index 370b19c..63db58b 100644 --- a/Basis/AtomicShellPair.ml +++ b/Basis/AtomicShellPair.ml @@ -28,12 +28,14 @@ let make ?(cutoff=Constants.epsilon) atomic_shell_a atomic_shell_b = let contracted_shell_pairs = List.map (fun s_a -> List.map (fun s_b -> - Csp.make ~cutoff s_a s_b + if Cs.index s_b <= Cs.index s_a then + Csp.make ~cutoff s_a s_b + else + None ) l_b ) l_a |> List.concat - |> List.filter (function None -> false | _ -> true) - |> List.map (function None -> assert false | Some x -> x) + |> list_some in match contracted_shell_pairs with | [] -> None diff --git a/Basis/AtomicShellPairCouple.ml b/Basis/AtomicShellPairCouple.ml index 6d90bae..e43078a 100644 --- a/Basis/AtomicShellPairCouple.ml +++ b/Basis/AtomicShellPairCouple.ml @@ -1,3 +1,5 @@ +open Util + type t = { atomic_shell_pair_p: AtomicShellPair.t ; @@ -32,8 +34,7 @@ let make ?(cutoff=Constants.epsilon) atomic_shell_pair_p atomic_shell_pair_q = ) (Asp.contracted_shell_pairs atomic_shell_pair_q) ) (Asp.contracted_shell_pairs atomic_shell_pair_p) |> List.concat - |> List.filter (function None -> false | _ -> true) - |> List.map (function None -> assert false | Some x -> x) + |> list_some in match contracted_shell_pair_couples with | [] -> None diff --git a/Basis/ContractedShell.ml b/Basis/ContractedShell.ml index 043e983..f8de68d 100644 --- a/Basis/ContractedShell.ml +++ b/Basis/ContractedShell.ml @@ -75,6 +75,8 @@ let size_of_shell x = Array.length x.norm_coef_scale let primitives x = x.prim +let zkey_array x = Ps.zkey_array x.prim.(0) + (** {2 Printers} *) diff --git a/Basis/ContractedShell.mli b/Basis/ContractedShell.mli index 865a4d3..4f87160 100644 --- a/Basis/ContractedShell.mli +++ b/Basis/ContractedShell.mli @@ -65,6 +65,9 @@ val norm_scales : t -> float array val size_of_shell : t -> int (** Number of contracted functions in the shell: length of {!norm_coef_scale}. *) +val zkey_array : t -> Zkey.t array +(** Returns the array of Zkeys associated with the contracted shell. *) + (** {2 Printers} *) diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index b92732b..e9dbd4b 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -120,11 +120,17 @@ let cmp a b = of contracted shells. *) let of_contracted_shell_array ?(cutoff=Constants.epsilon) basis = - Array.mapi (fun i shell_a -> - Array.mapi (fun j shell_b -> - make ~cutoff shell_a shell_b) - (Array.sub basis 0 (i+1)) - ) basis + let rec loop accu = function + | [] -> List.rev accu + | (s_a :: rest) as l -> + let new_accu = + (List.map (fun s_b -> make ~cutoff s_a s_b) l) :: accu + in loop new_accu rest + in + loop [] (Array.to_list basis) + |> List.concat + |> list_some + let equivalent x y = diff --git a/Basis/ContractedShellPair.mli b/Basis/ContractedShellPair.mli index f2c926d..30c0d91 100644 --- a/Basis/ContractedShellPair.mli +++ b/Basis/ContractedShellPair.mli @@ -22,18 +22,16 @@ val make : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option The contracted shell pair contains the only pairs of primitives for which the norm is greater than [cutoff]. - If all the primitive shell pairs are not significant, the function returns - [None]. + The function returns [None] if all the primitive shell pairs are not + significant, or if the index of the 1st primitive is smaller than the index + of the second primitive. + *) -val of_contracted_shell_array : ?cutoff:float -> ContractedShell.t array -> t option array array -(** Creates all possible contracted shell pairs from a list of contracted shells. - If a shell pair is not significant, sets the value to [None]: - -{[ - (of_contracted_shell_array p).(i).(j) = create p.(i) p.(j) -]} +val of_contracted_shell_array : ?cutoff:float -> ContractedShell.t array -> t list +(** Creates all possible contracted shell pairs from an array of contracted shells. + Only significant shell pairs are kept. *) val shell_a : t -> ContractedShell.t diff --git a/Basis/ContractedShellPairCouple.ml b/Basis/ContractedShellPairCouple.ml index b3c1b5f..2b47302 100644 --- a/Basis/ContractedShellPairCouple.ml +++ b/Basis/ContractedShellPairCouple.ml @@ -1,3 +1,5 @@ +open Util + type t = { shell_pair_p: ContractedShellPair.t ; @@ -35,8 +37,7 @@ let make ?(cutoff=Constants.epsilon) shell_pair_p shell_pair_q = ) (Csp.coefs_and_shell_pairs shell_pair_q) ) (Csp.coefs_and_shell_pairs shell_pair_p) |> List.concat - |> List.filter (function None -> false | _ -> true) - |> List.map (function None -> assert false | Some x -> x) + |> list_some in match coefs_and_shell_pair_couples with | [] -> None diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 76ef0f4..4d58689 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -7,6 +7,8 @@ open Bigarray type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t module Am = AngularMomentum +module As = AtomicShell +module Asp = AtomicShellPair module Bs = Basis module Cs = ContractedShell module Csp = ContractedShellPair @@ -49,6 +51,15 @@ let contracted_class_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : f TwoElectronRRVectorized.contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q +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 + +(* +let contracted_class_atomic_shell_pairs_vec ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = + TwoElectronRRVectorized.contracted_class_atomic_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q +*) + + let cutoff2 = cutoff *. cutoff (* type n_cls = { n : int ; cls : Zkey.t array } @@ -78,6 +89,9 @@ let of_basis basis = let n = Bs.size basis and shell = Bs.contracted_shells basis + (*TODO + and atomic_shells = Bs.atomic_shells basis + *) in @@ -85,34 +99,34 @@ let of_basis basis = let shell_pairs = Csp.of_contracted_shell_array shell in + (*TODO + let atomic_shell_pairs = + Asp.of_atomic_shell_array ~cutoff atomic_shells + in + *) (* Pre-compute diagonal integrals for Schwartz *) let t0 = Unix.gettimeofday () in let schwartz = - Array.map (fun pair_array -> Array.map (function - | None -> (Zmap.create 0, 0.) - | Some pair -> - let cls = - contracted_class_shell_pairs pair pair - in - (cls, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) - ) pair_array ) shell_pairs + List.map (fun pair -> + let cls = + contracted_class_shell_pairs pair pair + (*TODO + contracted_class_atomic_shell_pairs pair pair + *) + in + (pair, cls, 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) in - let icount = ref 0 in - for i=0 to (Array.length shell) - 1 do - print_int (Cs.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); + Printf.printf "%d shell pairs computed in %f seconds\n" + (List.length schwartz) (Unix.gettimeofday () -. t0); + (* Group shell pairs by common pairs of atoms *) - + (* 4D data initialization *) let eri_array = @@ -125,81 +139,86 @@ let of_basis basis = let t0 = Unix.gettimeofday () in let inn = ref 0 and out = ref 0 in - for i=0 to (Array.length shell) - 1 do - print_int (Cs.index shell.(i)) ; print_newline (); - for j=0 to i do - let schwartz_p, schwartz_p_max = schwartz.(i).(j) in + (*TODO + for i=0 to (Array.length atomic_shells) - 1 do + *) + let ishell = ref 0 in + List.iter (fun (shell_p, schwartz_p, schwartz_p_max) -> + 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 + (*TODO + Asp.atomic_shell_pairs shell_p + *) + in + try - if (schwartz_p_max < cutoff) then raise NullIntegral; + List.iter (fun (shell_q, schwartz_q, schwartz_q_max) -> + let () = + if Cs.index (Csp.shell_a shell_q) > + Cs.index (Csp.shell_a shell_p) then + raise Exit + in - let shell_p = - match shell_pairs.(i).(j) with - | None -> raise NullIntegral - | Some x -> x - in - - let sp = - Csp.shell_pairs shell_p - in - - for k=0 to i do - for l=0 to k do - let schwartz_q, schwartz_q_max = schwartz.(k).(l) in try - if schwartz_p_max *. schwartz_q_max < cutoff2 then - raise NullIntegral; + if schwartz_p_max *. schwartz_q_max < cutoff2 then + raise NullIntegral; - let shell_q = - match shell_pairs.(k).(l) with - | None -> raise NullIntegral - | Some x -> x - in + let sq = + Csp.shell_pairs shell_q + (*TODO + Asp.atomic_shell_pairs shell_q + *) + in - let sq = - Csp.shell_pairs shell_q - in - - let swap = + let swap = Array.length sp > Array.length sq - in + in - (* Compute all the integrals of the class *) - let cls = - if swap then - if (Array.length sp) + (Array.length sq) < 4 then - contracted_class_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p + (* Compute all the integrals of the class *) + let cls = + if swap then + (*TODO + contracted_class_atomic_shell_pairs ~schwartz_p:schwartz_q ~schwartz_q:schwartz_p shell_q shell_p + *) + if (Array.length sp) + (Array.length sq) < 4 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 - 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) < 4 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 + in - (* Write the data in the output file *) - Array.iteri (fun i_c powers_i -> - let i_c = Cs.index shell.(i) + i_c + 1 in - let xi = to_powers powers_i in - Array.iteri (fun j_c powers_j -> - let j_c = Cs.index shell.(j) + j_c + 1 in - let xj = to_powers powers_j in - Array.iteri (fun k_c powers_k -> - let k_c = Cs.index shell.(k) + k_c + 1 in - let xk = to_powers powers_k in - Array.iteri (fun l_c powers_l -> - let l_c = Cs.index shell.(l) + l_c + 1 in - let xl = to_powers powers_l in - let key = - if swap then - Zkey.of_powers_twelve xk xl xi xj - else - Zkey.of_powers_twelve xi xj xk xl - in - let value = - Zmap.find cls key - in + (* Write the data in the output file *) + 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 = + if swap then + Zkey.of_powers_twelve xk xl xi xj + else + Zkey.of_powers_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; @@ -208,21 +227,19 @@ let of_basis basis = 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; - ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(l)))) - ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(k)))) - ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(j)))) - ) Am.(zkey_array (Singlet (Cs.ang_mom shell.(i)))) + if (abs_float value > cutoff) then + (inn := !inn + 1; + ) + else + out := !out + 1; + ) (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)) with NullIntegral -> () - done; - done; - with NullIntegral -> () - done; - done; + ) schwartz + with Exit -> () + ) schwartz; Printf.printf "In: %d Out:%d\n" !inn !out ; Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0); eri_array diff --git a/Basis/PrimitiveShell.ml b/Basis/PrimitiveShell.ml index 48f280b..d709149 100644 --- a/Basis/PrimitiveShell.ml +++ b/Basis/PrimitiveShell.ml @@ -84,3 +84,5 @@ let norm_scales x = Lazy.force x.norm_scales let size_of_shell x = Array.length (norm_scales x) +let zkey_array x = Am.(zkey_array (Singlet (x.ang_mom))) + diff --git a/Basis/PrimitiveShell.mli b/Basis/PrimitiveShell.mli index b855c80..c0437a7 100644 --- a/Basis/PrimitiveShell.mli +++ b/Basis/PrimitiveShell.mli @@ -61,3 +61,6 @@ val norm_scales : t -> float array val size_of_shell : t -> int (** Number of functions in the shell. *) +val zkey_array : t -> Zkey.t array +(** Returns the array of Zkeys associated with the primitive shell. *) + diff --git a/Basis/PrimitiveShellPair.ml b/Basis/PrimitiveShellPair.ml index d575ceb..3646e1d 100644 --- a/Basis/PrimitiveShellPair.ml +++ b/Basis/PrimitiveShellPair.ml @@ -65,57 +65,57 @@ let create_make_of p_a p_b = function p_a -> - let norm_coef_a = - Ps.normalization p_a + let norm_coef_a = + Ps.normalization p_a + in + + let alpha_a = + Co.( Ps.exponent p_a |. Ps.center p_a ) + in + + function p_b -> + + let normalization = + norm_coef_a *. Ps.normalization p_b in - let alpha_a = - Co.( Ps.exponent p_a |. Ps.center p_a ) + let exponent = + Ps.exponent p_a +. Ps.exponent p_b in - function p_b -> + let exponent_inv = 1. /. exponent in - let normalization = - norm_coef_a *. Ps.normalization p_b + let normalization = + let argexpo = + Ps.exponent p_a *. Ps.exponent p_b *. a_minus_b_sq *. exponent_inv + in + normalization *. (pi *. exponent_inv)**1.5 *. exp (-. argexpo) + in + + function cutoff -> + + if abs_float normalization > cutoff then ( + + let beta_b = + Co.( Ps.exponent p_b |. Ps.center p_b ) in - let exponent = - Ps.exponent p_a +. Ps.exponent p_b + let center = + Co.(exponent_inv |. (alpha_a |+ beta_b)) in - let exponent_inv = 1. /. exponent in - - let normalization = - let argexpo = - Ps.exponent p_a *. Ps.exponent p_b *. a_minus_b_sq *. exponent_inv - in - normalization *. (pi *. exponent_inv)**1.5 *. exp (-. argexpo) + let center_minus_a = + Co.(center |- Ps.center p_a) in - function cutoff -> + Some { + ang_mom ; + exponent ; exponent_inv ; center ; center_minus_a ; a_minus_b ; + a_minus_b_sq ; normalization ; norm_scales ; shell_a = p_a; + shell_b = p_b } - if abs_float normalization > cutoff then ( - - let beta_b = - Co.( Ps.exponent p_b |. Ps.center p_b ) - in - - let center = - Co.(exponent_inv |. (alpha_a |+ beta_b)) - in - - let center_minus_a = - Co.(center |- Ps.center p_a) - in - - Some { - ang_mom ; - exponent ; exponent_inv ; center ; center_minus_a ; a_minus_b ; - a_minus_b_sq ; normalization ; norm_scales ; shell_a = p_a; - shell_b = p_b } - - ) - else None + ) + else None let make p_a p_b = diff --git a/Basis/PrimitiveShellPair.mli b/Basis/PrimitiveShellPair.mli index 8250ec6..e3cfb40 100644 --- a/Basis/PrimitiveShellPair.mli +++ b/Basis/PrimitiveShellPair.mli @@ -48,8 +48,7 @@ val create_make_of : PrimitiveShell.t -> PrimitiveShell.t -> are pre-computed. The result is None if the normalization coefficient of the resulting - function is below the cutoff, given as a last argument. - + function is below the cutoff given as a last argument. *) val ang_mom : t -> AngularMomentum.t diff --git a/Utils/Util.ml b/Utils/Util.ml index 1ac05d9..6641098 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -152,6 +152,15 @@ let boys_function ~maxm t = end +(** {2 List functions} *) + +let list_some l = + List.filter (function None -> false | _ -> true) l + |> List.map (function Some x -> x | _ -> assert false) + + +(** {2 Linear algebra} *) + let array_sum a = Array.fold_left ( +. ) 0. a diff --git a/Utils/Util.mli b/Utils/Util.mli index 27c42a4..4d99554 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -57,6 +57,10 @@ val array_product : float array -> float (** Returns the product of all the elements of the array *) +(** {2 Extension of the List module} *) +val list_some : 'a option list -> 'a list + + (** {2 Linear algebra } *) val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec