diff --git a/Basis/AtomicShell.ml b/Basis/AtomicShell.ml index ecc599e..438caa8 100644 --- a/Basis/AtomicShell.ml +++ b/Basis/AtomicShell.ml @@ -5,7 +5,7 @@ type t = { expo : float array array; coef : float array array; center : Coordinate.t; - totAngMom : AngularMomentum.t; + ang_mom : AngularMomentum.t; norm_coef : float array array; norm_coef_scale : float array; index : int; @@ -32,10 +32,10 @@ let make ?(index=0) contr = if not (unique_center (Array.length contr - 1)) then invalid_arg "ContractedAtomicShell.make Coordinate.t differ"; - let totAngMom = Cs.totAngMom contr.(0) in + let ang_mom = Cs.ang_mom contr.(0) in let rec unique_angmom = function | 0 -> true - | i -> if Cs.totAngMom contr.(i) = totAngMom then unique_angmom (i-1) else false + | i -> if Cs.ang_mom contr.(i) = ang_mom then unique_angmom (i-1) else false in if not (unique_angmom (Array.length contr - 1)) then invalid_arg "ContractedShell.make: AngularMomentum.t differ"; @@ -45,7 +45,7 @@ let make ?(index=0) contr = in let norm_coef_scale = Cs.norm_scales contr.(0) in - { index ; expo ; coef ; center ; totAngMom ; norm_coef ; + { index ; expo ; coef ; center ; ang_mom ; norm_coef ; norm_coef_scale ; contr } @@ -59,7 +59,7 @@ let coefficients x = x.coef let center x = x.center -let totAngMom x = x.totAngMom +let ang_mom x = x.ang_mom let size x = Array.length x.contr @@ -83,7 +83,7 @@ let pp_debug ppf x = fprintf ppf "@[<2>expo =@ %a ;@]@ " pp_float_2darray_size x.expo; fprintf ppf "@[<2>coef =@ %a ;@]@ " pp_float_2darray_size x.coef; fprintf ppf "@[<2>center =@ %a ;@]@ " Co.pp_angstrom x.center; - fprintf ppf "@[<2>totAngMom =@ %a ;@]@ " Am.pp_string x.totAngMom; + fprintf ppf "@[<2>ang_mom =@ %a ;@]@ " Am.pp_string x.ang_mom; fprintf ppf "@[<2>norm_coef =@ %a ;@]@ " pp_float_2darray_size x.norm_coef; fprintf ppf "@[<2>norm_coef_scale =@ %a ;@]@ " pp_float_array_size x.norm_coef_scale; fprintf ppf "@[<2>index =@ %d ;@]@ " x.index; @@ -91,7 +91,7 @@ let pp_debug ppf x = let pp ppf s = fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+ (size_of_shell s)*(size s)); - fprintf ppf "@[%a@ %a@] @[" Am.pp_string s.totAngMom Co.pp s.center; + fprintf ppf "@[%a@ %a@] @[" Am.pp_string s.ang_mom Co.pp s.center; Array.iter2 (fun e_arr c_arr -> fprintf ppf "@["; Array.iter2 (fun e c -> fprintf ppf "@[%16.8e %16.8e@]@;" e c) e_arr c_arr; diff --git a/Basis/AtomicShell.mli b/Basis/AtomicShell.mli index 7bab9c5..0610906 100644 --- a/Basis/AtomicShell.mli +++ b/Basis/AtomicShell.mli @@ -42,7 +42,7 @@ val index : t -> int val center : t -> Coordinate.t (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %}. *) -val totAngMom : t -> AngularMomentum.t +val ang_mom : t -> AngularMomentum.t (** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *) val size : t -> int @@ -68,7 +68,7 @@ val normalizations : t -> float array array val norm_scales : t -> float array (** Scaling factors {% $f(n_x,n_y,n_z)$ %}, given in the same order as - [AngularMomentum.zkey_array totAngMom]. *) + [AngularMomentum.zkey_array ang_mom]. *) val size_of_shell : t -> int (** Number of contracted functions in the shell: length of {!norm_coef_scale}. *) diff --git a/Basis/AtomicShellPair.mli b/Basis/AtomicShellPair.mli index cb3fdcc..1c6cb88 100644 --- a/Basis/AtomicShellPair.mli +++ b/Basis/AtomicShellPair.mli @@ -43,7 +43,7 @@ val norm_sq : t -> float val norm_scales : t -> float array (* norm_coef.(i) / norm_coef.(0) *) -val totAngMom : t -> AngularMomentum.t +val ang_mom : t -> AngularMomentum.t (* Total angular Momentum *) val monocentric : t -> bool diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 83058dc..6da5954 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -17,10 +17,10 @@ let of_nuclei_and_general_basis nucl bas = let contracted_shells = Array.map (fun (e, center) -> List.assoc e bas - |> Array.map (fun (totAngMom, shell) -> + |> Array.map (fun (ang_mom, shell) -> let lc = Array.map (fun Gb.{exponent ; coefficient} -> - coefficient, Ps.make totAngMom center exponent) shell + coefficient, Ps.make ang_mom center exponent) shell in let result = Cs.make ~index:!index_ lc in index_ := !index_ + Cs.size_of_shell result; @@ -32,16 +32,16 @@ let of_nuclei_and_general_basis nucl bas = in let atomic_shells = lazy( let uniq_center_angmom = - Array.map (fun x -> Cs.center x, Cs.totAngMom x) contracted_shells + Array.map (fun x -> Cs.center x, Cs.ang_mom x) contracted_shells |> Array.to_list |> List.sort_uniq compare in let csl = Array.to_list contracted_shells in - List.map (fun (center, totAngMom) -> + List.map (fun (center, ang_mom) -> let a = - List.filter (fun x -> Cs.center x = center && Cs.totAngMom x = totAngMom) csl + List.filter (fun x -> Cs.center x = center && Cs.ang_mom x = ang_mom) csl |> Array.of_list in As.make ~index:(Cs.index a.(0)) a diff --git a/Basis/ContractedShell.ml b/Basis/ContractedShell.ml index c648959..043e983 100644 --- a/Basis/ContractedShell.ml +++ b/Basis/ContractedShell.ml @@ -5,7 +5,7 @@ type t = { expo : float array; coef : float array; center : Coordinate.t; - totAngMom : AngularMomentum.t; + ang_mom : AngularMomentum.t; norm_coef : float array; norm_coef_scale : float array; index : int; @@ -32,10 +32,10 @@ let make ?(index=0) lc = if not (unique_center (Array.length prim - 1)) then invalid_arg "ContractedShell.make Coordinate.t differ"; - let totAngMom = Ps.totAngMom prim.(0) in + let ang_mom = Ps.ang_mom prim.(0) in let rec unique_angmom = function | 0 -> true - | i -> if Ps.totAngMom prim.(i) = totAngMom then unique_angmom (i-1) else false + | i -> if Ps.ang_mom prim.(i) = ang_mom then unique_angmom (i-1) else false in if not (unique_angmom (Array.length prim - 1)) then invalid_arg "ContractedShell.make: AngularMomentum.t differ"; @@ -47,7 +47,7 @@ let make ?(index=0) lc = in let norm_coef_scale = Ps.norm_scales prim.(0) in - { index ; expo ; coef ; center ; totAngMom ; norm_coef ; + { index ; expo ; coef ; center ; ang_mom ; norm_coef ; norm_coef_scale ; prim } @@ -61,7 +61,7 @@ let coefficients x = x.coef let center x = x.center -let totAngMom x = x.totAngMom +let ang_mom x = x.ang_mom let size x = Array.length x.prim @@ -85,18 +85,18 @@ let pp_debug ppf x = fprintf ppf "@[<2>expo =@ %a ;@]@ " pp_float_array_size x.expo; fprintf ppf "@[<2>coef =@ %a ;@]@ " pp_float_array_size x.coef; fprintf ppf "@[<2>center =@ %a ;@]@ " Co.pp_angstrom x.center; - fprintf ppf "@[<2>totAngMom =@ %a ;@]@ " Am.pp_string x.totAngMom; + fprintf ppf "@[<2>ang_mom =@ %a ;@]@ " Am.pp_string x.ang_mom; fprintf ppf "@[<2>norm_coef =@ %a ;@]@ " pp_float_array_size x.norm_coef; fprintf ppf "@[<2>norm_coef_scale =@ %a ;@]@ " pp_float_array_size x.norm_coef_scale; fprintf ppf "@[<2>index =@ %d ;@]@ " x.index; fprintf ppf "}@,@]" let pp ppf s = - (match s.totAngMom with + (match s.ang_mom with | Am.S -> fprintf ppf "@[%3d@] " (s.index+1) | _ -> fprintf ppf "@[%3d-%-3d@]" (s.index+1) (s.index+(Array.length s.norm_coef_scale)) ); - fprintf ppf "@[%a@ %a@]@[" Am.pp_string s.totAngMom Co.pp s.center; + fprintf ppf "@[%a@ %a@]@[" Am.pp_string s.ang_mom Co.pp s.center; Array.iter2 (fun e c -> fprintf ppf "@[%16.8e %16.8e@]@;" e c) s.expo s.coef; fprintf ppf "@]" diff --git a/Basis/ContractedShell.mli b/Basis/ContractedShell.mli index af0c1d1..865a4d3 100644 --- a/Basis/ContractedShell.mli +++ b/Basis/ContractedShell.mli @@ -40,7 +40,7 @@ val index : t -> int val center : t -> Coordinate.t (** Coordinate of the center {% $\mathbf{A} = (X_A,Y_A,Z_A)$ %}. *) -val totAngMom : t -> AngularMomentum.t +val ang_mom : t -> AngularMomentum.t (** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *) val size : t -> int @@ -60,7 +60,7 @@ val normalizations : t -> float array val norm_scales : t -> float array (** Scaling factors {% $f(n_x,n_y,n_z)$ %}, given in the same order as - [AngularMomentum.zkey_array totAngMom]. *) + [AngularMomentum.zkey_array ang_mom]. *) val size_of_shell : t -> int (** Number of contracted functions in the shell: length of {!norm_coef_scale}. *) diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index ece2116..3fde3ef 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -47,6 +47,7 @@ let make ?(cutoff=1.e-32) s_a s_b = | coefs_and_shell_pairs -> Some { shell_a = s_a ; shell_b = s_b ; coefs_and_shell_pairs } + let shell_a x = x.shell_a let shell_b x = x.shell_b let coefs_and_shell_pairs x = x.coefs_and_shell_pairs @@ -63,20 +64,20 @@ let exponents_inv x = List.map (fun (_,sp) -> Psp.exponent_inv sp) x.coefs_and_shell_pairs |> Array.of_list -let center_ab x = +let a_minus_b x = match x.coefs_and_shell_pairs with | [] -> assert false | (_,sp)::_ -> Psp.a_minus_b sp -let norm_sq x = +let a_minus_b_sq x = match x.coefs_and_shell_pairs with | [] -> assert false | (_,sp)::_ -> Psp.a_minus_b_sq sp -let totAngMom x = +let ang_mom x = match x.coefs_and_shell_pairs with | [] -> assert false - | (_,sp)::_ -> Psp.totAngMom sp + | (_,sp)::_ -> Psp.ang_mom sp let norm_scales x = match x.coefs_and_shell_pairs with @@ -89,6 +90,7 @@ let monocentric x = | (_,sp)::_ -> Psp.monocentric sp + (** Returns an integer characteristic of a contracted shell pair *) let hash a = Array.map Hashtbl.hash a @@ -151,22 +153,9 @@ let unique sp = aux [] sp -(** A map from a shell pair hash to the list of indices in the array - of shell pairs. -*) -let indices sp = - let map = - Hashtbl.create 129 - in - Array.iteri (fun i s -> - Array.iteri (fun j shell_p -> - let key = - hash shell_p - in - Hashtbl.add map key (i,j); ) s - ) sp; - map - - +let zkey_array x = + Am.zkey_array (Am.Doublet + Cs.(ang_mom x.shell_a, ang_mom x.shell_b) + ) diff --git a/Basis/ContractedShellPair.mli b/Basis/ContractedShellPair.mli index 2d2e108..faccce9 100644 --- a/Basis/ContractedShellPair.mli +++ b/Basis/ContractedShellPair.mli @@ -26,6 +26,7 @@ val make : ?cutoff:float -> ContractedShell.t -> ContractedShell.t -> t option [None]. *) + val of_contracted_shell_array : 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]: @@ -60,21 +61,25 @@ val coefficients : t -> float array val exponents_inv : t -> float array -val center_ab : t -> Coordinate.t +val a_minus_b : t -> Coordinate.t (* A-B *) -val norm_sq : t -> float +val a_minus_b_sq : t -> float (* |A-B|^2 *) val norm_scales : t -> float array (* normalizations.(i) / normalizations.(0) *) -val totAngMom : t -> AngularMomentum.t +val ang_mom : t -> AngularMomentum.t (* Total angular Momentum *) val monocentric : t -> bool (** If true, the two contracted shells have the same center. *) +val zkey_array : t -> Zkey.t array +(** Returns the array of Zkeys associated with the contracted shell pair. *) + + (* val hash : t -> int array val cmp : t -> t -> int diff --git a/Basis/ContractedShellPairCouple.ml b/Basis/ContractedShellPairCouple.ml new file mode 100644 index 0000000..3c6a7ce --- /dev/null +++ b/Basis/ContractedShellPairCouple.ml @@ -0,0 +1,75 @@ +type t = +{ + shell_pair_p: ContractedShellPair.t ; + shell_pair_q: ContractedShellPair.t ; + shell_a : ContractedShell.t ; + shell_b : ContractedShell.t ; + shell_c : ContractedShell.t ; + shell_d : ContractedShell.t ; + ang_mom : AngularMomentum.t ; + coefs_and_shell_pair_couples : (float * PrimitiveShellPairCouple.t) list ; +} + +module Am = AngularMomentum +module Co = Coordinate +module Cs = ContractedShell +module Csp = ContractedShellPair +module Pspc = PrimitiveShellPairCouple + +let make ?(cutoff=Constants.epsilon) shell_pair_p shell_pair_q = + let ang_mom = + Am.(Csp.ang_mom shell_pair_p + Csp.ang_mom shell_pair_q) + in + let shell_a = Csp.shell_a shell_pair_p + and shell_b = Csp.shell_b shell_pair_p + and shell_c = Csp.shell_a shell_pair_q + and shell_d = Csp.shell_b shell_pair_q + in + let cutoff = 1.e-3 *. cutoff in + let coefs_and_shell_pair_couples = + List.map (fun (c_ab, sp_ab) -> + List.map (fun (c_cd, sp_cd) -> + let coef_prod = c_ab *. c_cd in + if abs_float coef_prod < cutoff then None + else Some (coef_prod, Pspc.make sp_ab sp_cd) + ) (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) + in + match coefs_and_shell_pair_couples with + | [] -> None + | _ -> Some { shell_pair_p ; shell_pair_q ; ang_mom ; + shell_a ; shell_b ; shell_c ; shell_d ; + coefs_and_shell_pair_couples ; + } + +let coefs_and_shell_pair_couples t = t.coefs_and_shell_pair_couples + +let monocentric t = + Csp.monocentric t.shell_pair_p && Csp.monocentric t.shell_pair_q && + Csp.a_minus_b t.shell_pair_p = Csp.a_minus_b t.shell_pair_q + + +let ang_mom t = t.ang_mom + +let shell_pair_p t = t.shell_pair_p +let shell_pair_q t = t.shell_pair_q + +let shell_a t = t.shell_a +let shell_b t = t.shell_b +let shell_c t = t.shell_c +let shell_d t = t.shell_d + + +let zkey_array t = + match t.coefs_and_shell_pair_couples with + | (_,f)::_ -> Pspc.zkey_array f + | _ -> invalid_arg "ContractedShellPairCouple.zkey_array" + +let norm_scales t = + match t.coefs_and_shell_pair_couples with + | (_,f)::_ -> Pspc.norm_scales f + | _ -> invalid_arg "ContractedShellPairCouple.norm_scales" + diff --git a/Basis/ContractedShellPairCouple.mli b/Basis/ContractedShellPairCouple.mli new file mode 100644 index 0000000..4de4016 --- /dev/null +++ b/Basis/ContractedShellPairCouple.mli @@ -0,0 +1,61 @@ +(** Data structure describing a couple of contracted shells pairs. + +A contracted shell pair couple is the cartesian product between two sets of functions, one +set over electron one and one set over electron two. Both sets are contracted shell +pairs. + +These are usually called {e shell quartets} in the literature, but we prefer to use +{e pair} for two functions with the same electron, and {e couple} for two functions +acting on different electrons, since they will be coupled by a two-electron operator. + + +*) + +type t + + +val make : ?cutoff:float -> ContractedShellPair.t -> ContractedShellPair.t -> t option +(** Creates a contracted shell pair couple using two contracted shell pairs. +*) + +val ang_mom : t -> AngularMomentum.t +(** Total angular momentum of the shell pair couple: sum of the angular momenta of + all the shells. *) + + +val shell_a : t -> ContractedShell.t +(** Returns the first contracted shell of the first shell pair. *) + +val shell_b : t -> ContractedShell.t +(** Returns the second contracted shell of the first shell pair. *) + +val shell_c : t -> ContractedShell.t +(** Returns the first contracted shell of the second shell pair. *) + +val shell_d : t -> ContractedShell.t +(** Returns the second contracted shell of the second shell pair. *) + + +val shell_pair_p : t -> ContractedShellPair.t +(** Returns the first contracted shell pair that was used to build the shell pair. *) + +val shell_pair_q : t -> ContractedShellPair.t +(** Returns the second contracted shell pair that was used to build the shell pair. *) + +val monocentric : t -> bool +(** True if all four contracted shells have the same center. *) + +val coefs_and_shell_pair_couples : t -> (float * PrimitiveShellPairCouple.t) list +(** Returns the list of significant shell pair couples. *) + +val zkey_array : t -> Zkey.t array +(** Returns the array of {!Zkey.t} relative to the four shells of the + shell pair couple. +*) + +val norm_scales : t -> float array +(** Scaling factors of normalization coefficients inside the shell. The + ordering is the same as {!zkey_array}. +*) + + diff --git a/Basis/ERI.ml b/Basis/ERI.ml index ad05f22..bacbb65 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -221,10 +221,10 @@ let of_basis basis = ) else out := !out + 1; - ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(l)))) - ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(k)))) - ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(j)))) - ) Am.(zkey_array (Singlet (Cs.totAngMom shell.(i)))) + ) 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)))) with NullIntegral -> () done; done; @@ -280,7 +280,7 @@ let to_file ~filename basis = let oc = open_out filename in let zkey = Array.map (fun b -> let result = - Angular_momentum.(zkey_array (Kind_1 b.totAngMom)) + Angular_momentum.(zkey_array (Kind_1 b.ang_mom)) in { n=Array.length result ; cls=result } ) basis @@ -378,7 +378,7 @@ let to_file ~filename basis = let xto_file ~filename basis = let zkey = Array.map (fun b -> let result = - Angular_momentum.(zkey_array (Kind_1 b.totAngMom)) + Angular_momentum.(zkey_array (Kind_1 b.ang_mom)) in { n=Array.length result ; cls=result } ) basis diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index f701e9f..c5186ea 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -28,9 +28,7 @@ let contracted_class shell_a shell_b : float Zmap.t = | Some shell_p -> begin (* Pre-computation of integral class indices *) - let class_indices = - Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b)) - in + let class_indices = Csp.zkey_array shell_p in let contracted_class = Array.make (Array.length class_indices) 0. @@ -39,8 +37,8 @@ let contracted_class shell_a shell_b : float Zmap.t = (* Compute all integrals in the shell for each pair of significant shell pairs *) let sp = Csp.shell_pairs shell_p in - let center_ab = - Csp.center_ab shell_p + let a_minus_b = + Csp.a_minus_b shell_p in let norm_coef_scales = Csp.norm_scales shell_p @@ -78,7 +76,7 @@ let contracted_class shell_a shell_b : float Zmap.t = let xyz = xyz_of_int k in Overlap_primitives.hvrr (a, b) expo_inv - (Co.get xyz center_ab, + (Co.get xyz a_minus_b, Co.get xyz center_pa) in let f k = @@ -154,8 +152,8 @@ let of_basis basis = in result.{i_c,j_c} <- value; result.{j_c,i_c} <- value; - ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i)))) - ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) done; done; Mat.detri result; diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml index c49390a..9b555e7 100644 --- a/Basis/NucInt.ml +++ b/Basis/NucInt.ml @@ -87,8 +87,8 @@ let of_basis_nuclei basis nuclei = in eni_array.{j_c,i_c} <- value; eni_array.{i_c,j_c} <- value; - ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j)))) - ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i)))) done; done; Mat.detri eni_array; diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index e1103c1..96bb9cb 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -101,15 +101,10 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let shell_a = Csp.shell_a shell_p and shell_b = Csp.shell_b shell_p in - let maxm = - Am.(Cs.totAngMom shell_a + Cs.totAngMom shell_b) - |> Am.to_int - in + let maxm = Am.to_int (Csp.ang_mom shell_p) in (* Pre-computation of integral class indices *) - let class_indices = - Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b)) - in + let class_indices = Csp.zkey_array shell_p in let contracted_class = Array.make (Array.length class_indices) 0.; @@ -119,7 +114,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let norm_scales_p = Csp.norm_scales shell_p in - let center_ab = Csp.center_ab shell_p in + let center_ab = Csp.a_minus_b shell_p in List.iter (fun (coef_prod, psp) -> try @@ -146,7 +141,7 @@ let contracted_class_shell_pair ~zero_m shell_p geometry : float Zmap.t = let zero_m_array = zero_m ~maxm ~expo_pq_inv ~norm_pq_sq in - match Cs.(totAngMom shell_a, totAngMom shell_b) with + match Cs.(ang_mom shell_a, ang_mom shell_b) with | Am.(S,S) -> let integral = zero_m_array.(0) in contracted_class.(0) <- contracted_class.(0) -. coef_prod *. integral *. charge diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index 5447a1b..2e98a8d 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -28,16 +28,14 @@ let contracted_class shell_a shell_b : float Zmap.t = begin (* Pre-computation of integral class indices *) - let class_indices = - Am.zkey_array (Am.Doublet Cs.(totAngMom shell_a, totAngMom shell_b)) - in + let class_indices = Csp.zkey_array shell_p in let contracted_class = Array.make (Array.length class_indices) 0. in - let center_ab = - Csp.center_ab shell_p + let a_minus_b = + Csp.a_minus_b shell_p in let norm_coef_scales = Csp.norm_scales shell_p @@ -66,7 +64,7 @@ let contracted_class shell_a shell_b : float Zmap.t = let xyz = xyz_of_int k in Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB) expo_inv - (Co.get xyz center_ab, + (Co.get xyz a_minus_b, Co.get xyz center_pa) in let norm = norm_coef_scales.(i) in @@ -120,8 +118,8 @@ let of_basis basis = in result.{i_c,j_c} <- value; result.{j_c,i_c} <- value; - ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(i)))) - ) (Am.zkey_array (Singlet (Cs.totAngMom shell.(j)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) done; done; Mat.detri result; diff --git a/Basis/PrimitiveShell.ml b/Basis/PrimitiveShell.ml index 00d8d39..48f280b 100644 --- a/Basis/PrimitiveShell.ml +++ b/Basis/PrimitiveShell.ml @@ -7,15 +7,15 @@ type t = { normalization : float; norm_scales : float array lazy_t; center : Coordinate.t; - totAngMom : AngularMomentum.t; + ang_mom : AngularMomentum.t; } module Am = AngularMomentum -let compute_norm_coef alpha totAngMom = +let compute_norm_coef alpha ang_mom = let atot = - Am.to_int totAngMom + Am.to_int ang_mom in let factor int_array = let dfa = Array.map (fun j -> @@ -37,15 +37,15 @@ let compute_norm_coef alpha totAngMom = in f -let make totAngMom center exponent = +let make ang_mom center exponent = let norm_coef_func = - compute_norm_coef exponent totAngMom + compute_norm_coef exponent ang_mom in let norm = - 1. /. norm_coef_func [| Am.to_int totAngMom ; 0 ; 0 |] + 1. /. norm_coef_func [| Am.to_int ang_mom ; 0 ; 0 |] in let powers = - Am.zkey_array (Am.Singlet totAngMom) + Am.zkey_array (Am.Singlet ang_mom) in let norm_scales = lazy ( Array.map (fun a -> @@ -53,12 +53,12 @@ let make totAngMom center exponent = ) powers ) in let normalization = 1. /. norm in - { exponent ; normalization ; norm_scales ; center ; totAngMom } + { exponent ; normalization ; norm_scales ; center ; ang_mom } let to_string s = let coord = s.center in - Printf.sprintf "%1s %8.3f %8.3f %8.3f %16.8e" (Am.to_string s.totAngMom) + Printf.sprintf "%1s %8.3f %8.3f %8.3f %16.8e" (Am.to_string s.ang_mom) (get X coord) (get Y coord) (get Z coord) s.exponent @@ -74,7 +74,7 @@ let exponent x = x.exponent let center x = x.center -let totAngMom x = x.totAngMom +let ang_mom x = x.ang_mom let norm x = 1. /. x.normalization diff --git a/Basis/PrimitiveShell.mli b/Basis/PrimitiveShell.mli index 5a400fa..b855c80 100644 --- a/Basis/PrimitiveShell.mli +++ b/Basis/PrimitiveShell.mli @@ -33,7 +33,7 @@ val exponent : t -> float val center : t -> Coordinate.t (** Coordinate {% $\mathbf{A}$ %}.of the center. *) -val totAngMom : t -> AngularMomentum.t +val ang_mom : t -> AngularMomentum.t (** Total angular momentum : {% $l = n_x + n_y + n_z$ %}. *) val norm : t -> float @@ -54,7 +54,7 @@ val norm_scales : t -> float array (** Scaling factors {% $f(n_x,n_y,n_z)$ %} adjusting the normalization coefficient for the powers of {% $x,y,z$ %}. The normalization coefficients of the functions of the shell are given by {% $\mathcal{N}\times f$ %}. They are - given in the same order as [AngularMomentum.zkey_array totAngMom]: + given in the same order as [AngularMomentum.zkey_array ang_mom]: {% \\[ f(n_x,n_y,n_z) = \frac{|| g_{l,0,0}(\mathbf{r}) ||}{|| g_{n_x,n_y,n_z}(\mathbf{r}) ||} \\] %} *) diff --git a/Basis/PrimitiveShellPair.ml b/Basis/PrimitiveShellPair.ml index 3882938..d575ceb 100644 --- a/Basis/PrimitiveShellPair.ml +++ b/Basis/PrimitiveShellPair.ml @@ -12,7 +12,7 @@ type t = { norm_scales : float array lazy_t; normalization : float; (* norm_coef_a * norm_coef_b * g, with g = (pi/(alpha+beta))^(3/2) exp (-|A-B|^2 * alpha*beta/(alpha+beta)) *) - totAngMom : AngularMomentum.t; + ang_mom : AngularMomentum.t; shell_a : PrimitiveShell.t; shell_b : PrimitiveShell.t; } @@ -30,7 +30,7 @@ let hash a = let equivalent a b = a.exponent = b.exponent && - a.totAngMom = b.totAngMom && + a.ang_mom = b.ang_mom && a.normalization = b.normalization && a.center = b.center && a.center_minus_a = b.center_minus_a && @@ -59,8 +59,8 @@ let create_make_of p_a p_b = |> Array.concat ) in - let totAngMom = - Am.( Ps.totAngMom p_a + Ps.totAngMom p_b ) + let ang_mom = + Am.( Ps.ang_mom p_a + Ps.ang_mom p_b ) in function p_a -> @@ -109,7 +109,7 @@ let create_make_of p_a p_b = in Some { - totAngMom ; + 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 } @@ -136,7 +136,7 @@ let monocentric x = Ps.center x.shell_a = Ps.center x.shell_b -let totAngMom x = x.totAngMom +let ang_mom x = x.ang_mom let a_minus_b x = x.a_minus_b @@ -154,3 +154,11 @@ let shell_a x = x.shell_a let shell_b x = x.shell_b + +let zkey_array x = + Am.zkey_array (Am.Doublet + Ps.(ang_mom x.shell_a, ang_mom x.shell_b) + ) + + + diff --git a/Basis/PrimitiveShellPair.mli b/Basis/PrimitiveShellPair.mli index 42c3e67..8250ec6 100644 --- a/Basis/PrimitiveShellPair.mli +++ b/Basis/PrimitiveShellPair.mli @@ -1,7 +1,7 @@ (** Data structure describing a pair of primitive shells. A primitive shell pair is the cartesian product between two sets of functions, each -set containing all the functions of a primitive shell. +set containing all the functions of a primitive shell. These are one-electron functions. {% \\[ \left\\{ p_{k_x,k_y,k_z}(\mathbf{r}) \right\\} = @@ -24,7 +24,7 @@ where {!a_minus_b}, {!a_minus_b_sq} and {!norm_coef_scale} depend only on the -centering of the two shells, and {!totAngMom} only depends on the angular +centering of the two shells, and {!ang_mom} only depends on the angular momenta of the two shells. Hence, these quantities need to be computed only once when a {!ContractedShellPair.t} is built. Hence, there is the {!create_make_of} function which creates a [make] function which is suitable @@ -52,7 +52,7 @@ val create_make_of : PrimitiveShell.t -> PrimitiveShell.t -> *) -val totAngMom : t -> AngularMomentum.t +val ang_mom : t -> AngularMomentum.t (** Total angular momentum of the shell pair: sum of the angular momenta of the shells. *) @@ -104,3 +104,6 @@ val hash : t -> int val cmp : t -> t -> int (** Arbitray comparison function for sorting. *) +val zkey_array : t -> Zkey.t array +(** Returns the array of Zkeys associated with the shell pair. *) + diff --git a/Basis/PrimitiveShellPairCouple.ml b/Basis/PrimitiveShellPairCouple.ml new file mode 100644 index 0000000..70dc88c --- /dev/null +++ b/Basis/PrimitiveShellPairCouple.ml @@ -0,0 +1,71 @@ +exception NullContribution + +type t = +{ + shell_pair_p: PrimitiveShellPair.t ; + shell_pair_q: PrimitiveShellPair.t ; + shell_a : PrimitiveShell.t ; + shell_b : PrimitiveShell.t ; + shell_c : PrimitiveShell.t ; + shell_d : PrimitiveShell.t ; + ang_mom : AngularMomentum.t ; + zkey_array : Zkey.t array lazy_t; +} + +module Am = AngularMomentum +module Co = Coordinate +module Ps = PrimitiveShell +module Psp = PrimitiveShellPair + +let make shell_pair_p shell_pair_q = + let ang_mom = + Am.(Psp.ang_mom shell_pair_p + Psp.ang_mom shell_pair_q) + in + let shell_a = Psp.shell_a shell_pair_p + and shell_b = Psp.shell_b shell_pair_p + and shell_c = Psp.shell_a shell_pair_q + and shell_d = Psp.shell_b shell_pair_q + in + let zkey_array = lazy ( + Am.zkey_array (Am.Quartet + Ps.(ang_mom shell_a, ang_mom shell_b, + ang_mom shell_c, ang_mom shell_d) + ) + ) + in + { shell_pair_p ; shell_pair_q ; ang_mom ; + shell_a ; shell_b ; shell_c ; shell_d ; + zkey_array ; + } + + + +let monocentric t = + Psp.monocentric t.shell_pair_p && Psp.monocentric t.shell_pair_q && + Psp.center t.shell_pair_p = Psp.center t.shell_pair_q + + +let ang_mom t = t.ang_mom + +let shell_pair_p t = t.shell_pair_p +let shell_pair_q t = t.shell_pair_q + +let shell_a t = t.shell_a +let shell_b t = t.shell_b +let shell_c t = t.shell_c +let shell_d t = t.shell_d + +let p_minus_q t = + let p = Psp.center t.shell_pair_p + and q = Psp.center t.shell_pair_q + in Co.(p |- q) + +let zkey_array t = Lazy.force t.zkey_array + +let norm_scales t = + let norm_coef_scale_p_list = Array.to_list (Psp.norm_scales t.shell_pair_p) in + let norm_coef_scale_q = Psp.norm_scales t.shell_pair_q in + List.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) + norm_coef_scale_p_list + |> Array.concat + diff --git a/Basis/PrimitiveShellPairCouple.mli b/Basis/PrimitiveShellPairCouple.mli new file mode 100644 index 0000000..cce84f3 --- /dev/null +++ b/Basis/PrimitiveShellPairCouple.mli @@ -0,0 +1,62 @@ +(** Data structure describing a couple of primitive shells pairs. + +A primitive shell pair couple is the cartesian product between two sets of functions, one +set over electron one and one set over electron two. Both sets are primitive shell +pairs. + +These are usually called {e shell quartets} in the literature, but we prefer to use +{e pair} for two functions with the same electron, and {e couple} for two functions +acting on different electrons, since they will be coupled by a two-electron operator. + + +*) + +type t + + +val make : PrimitiveShellPair.t -> PrimitiveShellPair.t -> t +(** Creates a primitive shell pair couple using two primitive shell pairs. +*) + +val ang_mom : t -> AngularMomentum.t +(** Total angular momentum of the shell pair couple: sum of the angular momenta of + all the shells. *) + + +val shell_a : t -> PrimitiveShell.t +(** Returns the first primitive shell of the first shell pair. *) + +val shell_b : t -> PrimitiveShell.t +(** Returns the second primitive shell of the first shell pair. *) + +val shell_c : t -> PrimitiveShell.t +(** Returns the first primitive shell of the second shell pair. *) + +val shell_d : t -> PrimitiveShell.t +(** Returns the second primitive shell of the second shell pair. *) + + +val shell_pair_p : t -> PrimitiveShellPair.t +(** Returns the first primitive shell pair that was used to build the shell pair. *) + +val shell_pair_q : t -> PrimitiveShellPair.t +(** Returns the second primitive shell pair that was used to build the shell pair. *) + +val p_minus_q : t -> Coordinate.t +(** {% $\mathbf{P-Q}$ %}. *) + + +val monocentric : t -> bool +(** True if all four primitive shells have the same center. *) + +val zkey_array : t -> Zkey.t array +(** Returns the array of {!Zkey.t} relative to the four shells of the + shell pair couple. +*) + +val norm_scales : t -> float array +(** Scaling factors of normalization coefficients inside the shell. The + ordering is the same as {!zkey_array}. +*) + + diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index fdcbbb0..4d37f05 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -1,15 +1,18 @@ open Util -module Am = AngularMomentum -module Co = Coordinate -module Cs = ContractedShell -module Csp = ContractedShellPair -module Po = Powers -module Psp = PrimitiveShellPair -module Ps = PrimitiveShell +module Am = AngularMomentum +module Co = Coordinate +module Cs = ContractedShell +module Csp = ContractedShellPair +module Cspc = ContractedShellPairCouple +module Po = Powers +module Psp = PrimitiveShellPair +module Pspc = PrimitiveShellPairCouple +module Ps = PrimitiveShell let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff +let empty = Zmap.create 0 exception NullQuartet @@ -275,50 +278,39 @@ let rec hvrr_two_e let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - let shell_a = Csp.shell_a shell_p - and shell_b = Csp.shell_b shell_p - and shell_c = Csp.shell_a shell_q - and shell_d = Csp.shell_b shell_q - in - let maxm = Am.(Csp.totAngMom shell_p + Csp.totAngMom shell_q |> to_int) in + match Cspc.make ~cutoff shell_p shell_q with + | None -> empty + | Some shell_pair_couple -> - (* Pre-computation of integral class indices *) - let class_indices = - Am.zkey_array (Am.Quartet - Cs.(totAngMom shell_a, totAngMom shell_b, - totAngMom shell_c, totAngMom shell_d )) - in + let maxm = Am.to_int (Cspc.ang_mom shell_pair_couple) in - let contracted_class = - Array.make (Array.length class_indices) 0.; - in + (* Pre-computation of integral class indices *) + let class_indices = Cspc.zkey_array shell_pair_couple in - let monocentric = - Csp.monocentric shell_p && Csp.monocentric shell_q && - Cs.center (Csp.shell_a shell_p) = Cs.center (Csp.shell_a shell_q) - in - - (* Compute all integrals in the shell for each pair of significant shell pairs *) - - let norm_coef_scale_p_list = Array.to_list (Csp.norm_scales shell_p) in - let norm_coef_scale_q = Csp.norm_scales shell_q in - - let center_ab = Csp.center_ab shell_p in - List.iter (fun (c_ab, sp_ab) -> - - let expo_b = Ps.exponent (Psp.shell_b sp_ab) - and expo_inv_p = Psp.exponent_inv sp_ab - and center_pa = Psp.center_minus_a sp_ab + let contracted_class = + Array.make (Array.length class_indices) 0.; in - List.iter (fun (c_cd, sp_cd) -> + let monocentric = + Cspc.monocentric shell_pair_couple + in - let coef_prod = c_ab *. c_cd in + (* Compute all integrals in the shell for each pair of significant shell pairs *) - (** Screening on the product of coefficients *) - try - if (abs_float coef_prod) < 1.e-3 *. cutoff then - raise NullQuartet; + let center_ab = Csp.a_minus_b shell_p + and center_cd = Csp.a_minus_b shell_q + in + + let norm_scales = Cspc.norm_scales shell_pair_couple in + + List.iter (fun (coef_prod, spc) -> + + let sp_ab = Pspc.shell_pair_p spc + and sp_cd = Pspc.shell_pair_q spc + in + + let expo_inv_p = Psp.exponent_inv sp_ab + in let center_pq = Co.(Psp.center sp_ab |- Psp.center sp_cd) in let norm_pq_sq = Co.dot center_pq center_pq in @@ -330,23 +322,19 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in begin - match Cs.(totAngMom shell_a, totAngMom shell_b, - totAngMom shell_c, totAngMom shell_d) with - | Am.(S,S,S,S) -> - let integral = - zero_m_array.(0) - in + match Cspc.ang_mom shell_pair_couple with + | Am.S -> + let integral = zero_m_array.(0) in contracted_class.(0) <- contracted_class.(0) +. coef_prod *. integral | _ -> - let expo_d = Ps.exponent (Psp.shell_b sp_cd) in - let map_1d = Zmap.create (4*maxm) in - let map_2d = Zmap.create (Array.length class_indices) in - let center_cd = Psp.a_minus_b sp_cd in - let center_qc = Psp.center_minus_a sp_cd in - let norm_coef_scale = - List.map (fun v1 -> Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) - norm_coef_scale_p_list - |> Array.concat + let expo_b = Ps.exponent (Psp.shell_b sp_ab) + and expo_d = Ps.exponent (Psp.shell_b sp_cd) + and center_pa = Psp.center_minus_a sp_ab + in + let map_1d = Zmap.create (4*maxm) + and map_2d = Zmap.create (Array.length class_indices) + in + let center_qc = Psp.center_minus_a sp_cd in (* Compute the integral class from the primitive shell quartet *) @@ -363,7 +351,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q if ( ((1 land angMom_a.Po.x + angMom_b.Po.x + angMom_c.Po.x + angMom_d.Po.x)=1) || ((1 land angMom_a.Po.y + angMom_b.Po.y + angMom_c.Po.y + angMom_d.Po.y)=1) || ((1 land angMom_a.Po.z + angMom_b.Po.z + angMom_c.Po.z + angMom_d.Po.z)=1) - ) then + ) then raise NullQuartet end; @@ -393,7 +381,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q *) - let norm = norm_coef_scale.(i) in + let norm = norm_scales.(i) in let coef_prod = coef_prod *. norm in let integral = @@ -410,15 +398,13 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q with NullQuartet -> () ) end - with NullQuartet -> () - ) (Csp.coefs_and_shell_pairs shell_q) - ) (Csp.coefs_and_shell_pairs shell_p); + ) (Cspc.coefs_and_shell_pair_couples shell_pair_couple); - let result = - Zmap.create (Array.length contracted_class) - in - Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; - result + let result = + Zmap.create (Array.length contracted_class) + in + Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; + result @@ -431,6 +417,6 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = match shell_p, shell_q with | Some shell_p, Some shell_q -> contracted_class_shell_pairs ~zero_m shell_p shell_q - | _ -> Zmap.create 0 + | _ -> empty diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 3ba4031..7a33819 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -2,19 +2,21 @@ open Util open Lacaml.D open Bigarray -module Am = AngularMomentum -module Co = Coordinate -module Cs = ContractedShell -module Csp = ContractedShellPair -module Po = Powers -module Psp = PrimitiveShellPair -module Ps = PrimitiveShell +module Am = AngularMomentum +module Co = Coordinate +module Cs = ContractedShell +module Csp = ContractedShellPair +module Cspc = ContractedShellPairCouple +module Po = Powers +module Psp = PrimitiveShellPair +module Ps = PrimitiveShell exception NullQuartet exception Found let cutoff = Constants.integrals_cutoff let cutoff2 = cutoff *. cutoff +let empty = Zmap.create 0 let at_least_one_valid arr = try @@ -555,57 +557,10 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q : float Zmap.t = - let shell_a = Csp.shell_a shell_p - and shell_b = Csp.shell_b shell_p - and shell_c = Csp.shell_a shell_q - and shell_d = Csp.shell_b shell_q - in - let maxm = - Am.(Csp.totAngMom shell_p + Csp.totAngMom shell_q |> to_int) - in - - (* Pre-computation of integral class indices *) - let class_indices = - Am.zkey_array (Am.Quartet - Cs.(totAngMom shell_a, totAngMom shell_b, - totAngMom shell_c, totAngMom shell_d)) - in - - let contracted_class = - Array.make (Array.length class_indices) 0.; - in - - let monocentric = - Csp.monocentric shell_p && Csp.monocentric shell_q && - Cs.center (Csp.shell_a shell_p) = - Cs.center (Csp.shell_a shell_q) - in - - (** Screening on the product of coefficients *) - let coef_max_p = - List.fold_left (fun accu (x,_) -> - if (abs_float x) > accu then (abs_float x) else accu) - 0. (Csp.coefs_and_shell_pairs shell_p) - and coef_max_q = - List.fold_left (fun accu (x,_) -> - if (abs_float x) > accu then (abs_float x) else accu) - 0. (Csp.coefs_and_shell_pairs shell_q) - in - - let filter_p = - let cutoff_p = cutoff /. coef_max_p in - Csp.coefs_and_shell_pairs shell_p - |> List.filter (fun (ck,f) -> abs_float ck > cutoff_p) - and filter_q = - let cutoff_q = cutoff /. coef_max_q in - Csp.coefs_and_shell_pairs shell_q - |> List.filter (fun (ck,f) -> abs_float ck > cutoff_q) - in - - let sp = List.map snd filter_p |> Array.of_list - and sq = List.map snd filter_q |> Array.of_list - and cp = List.map fst filter_p |> Array.of_list - and cq = List.map fst filter_q |> Array.of_list + let sp = Csp.shell_pairs shell_p + and sq = Csp.shell_pairs shell_q + and cp = Csp.coefficients shell_p + and cq = Csp.coefficients shell_q in let np, nq = @@ -613,207 +568,219 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q Array.length sq in + try + match Cspc.make ~cutoff shell_p shell_q with + | None -> raise NullQuartet + | Some shell_pair_couple -> - (* Compute all integrals in the shell for each pair of significant shell pairs *) - - begin - match Cs.(totAngMom shell_a, totAngMom shell_b, - totAngMom shell_c, totAngMom shell_d) with - | Am.(S,S,S,S) -> - contracted_class.(0) <- - begin - try - let expo_inv_p = - Vec.init np (fun ab -> Psp.exponent_inv sp.(ab-1)) - and expo_inv_q = - Vec.init nq (fun cd -> Psp.exponent_inv sq.(cd-1)) - in - - let coef = - let result = Mat.make0 nq np in - Lacaml.D.ger (Vec.of_array @@ cq) (Vec.of_array @@ cp) result; - result - in - let zm_array = Mat.init_cols np nq (fun i j -> - try - if (abs_float coef.{j,i} ) < 1.e-3*.cutoff then - raise NullQuartet; - - let expo_pq_inv = - expo_inv_p.{i} +. expo_inv_q.{j} - in - - let center_pq = - Co.(Psp.center sp.(i-1) |- Psp.center sq.(j-1)) - in - let norm_pq_sq = - Co.dot center_pq center_pq - in - - let zero_m_array = - zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq - in - zero_m_array.(0) - with NullQuartet -> 0. - ) in - Mat.gemm_trace zm_array coef - with (Invalid_argument _) -> 0. - end - | _ -> - - let coef = - Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) ) + let shell_a = Cspc.shell_a shell_pair_couple + and shell_c = Cspc.shell_c shell_pair_couple in - let expo_inv_p = - Array.map (fun shell_ab -> Psp.exponent_inv shell_ab) sp - and expo_inv_q = - Array.map (fun shell_cd -> Psp.exponent_inv shell_cd) sq + let maxm = Am.to_int (Cspc.ang_mom shell_pair_couple) in + + + (* Pre-computation of integral class indices *) + let class_indices = Cspc.zkey_array shell_pair_couple in + + let contracted_class = + Array.make (Array.length class_indices) 0.; in - let expo_b = - Array.map (fun shell_ab -> Ps.exponent (Psp.shell_b shell_ab) ) sp - and expo_d = - Array.map (fun shell_cd -> Ps.exponent (Psp.shell_b shell_cd) ) sq + let monocentric = + Cspc.monocentric shell_pair_couple in - let norm_coef_scale_p = Csp.norm_scales shell_p in - let center_pq = - let result = - Array.init 3 (fun xyz -> - Array.init np (fun ab -> - let shell_ab = sp.(ab) in - Array.init nq (fun cd -> - let shell_cd = sq.(cd) + (* Compute all integrals in the shell for each pair of significant shell pairs *) + + begin + match Cspc.ang_mom shell_pair_couple with + | Am.S -> + contracted_class.(0) <- + begin + try + let expo_inv_p = + Vec.init np (fun ab -> Psp.exponent_inv sp.(ab-1)) + and expo_inv_q = + Vec.init nq (fun cd -> Psp.exponent_inv sq.(cd-1)) in - let cpq = - Co.(Psp.center shell_ab |- Psp.center shell_cd) + + let coef = + let result = Mat.make0 nq np in + Lacaml.D.ger (Vec.of_array @@ cq) (Vec.of_array @@ cp) result; + result in - match xyz with - | 0 -> Co.get X cpq; - | 1 -> Co.get Y cpq; - | _ -> Co.get Z cpq; - ) - ) - ) - in function - | Co.X -> result.(0) - | Co.Y -> result.(1) - | Co.Z -> result.(2) - in - let center_pa = - let result = - Array.init 3 (fun xyz -> - Array.init np (fun ab -> - let shell_ab = sp.(ab) in - let cpa = - Co.(Psp.center shell_ab |- Cs.center shell_a) - in - match xyz with - | 0 -> Co.(get X cpa); - | 1 -> Co.(get Y cpa); - | _ -> Co.(get Z cpa); - ) - ) - in function - | Co.X -> result.(0) - | Co.Y -> result.(1) - | Co.Z -> result.(2) - in - let center_qc = - let result = - Array.init 3 (fun xyz -> - Array.init nq (fun cd -> - let shell_cd = sq.(cd) in - let cqc = - Co.(Psp.center shell_cd |- Cs.center shell_c) - in - match xyz with - | 0 -> Co.(get X cqc); - | 1 -> Co.(get Y cqc); - | _ -> Co.(get Z cqc); - ) - ) - in function - | Co.X -> result.(0) - | Co.Y -> result.(1) - | Co.Z -> result.(2) - in - let zero_m_array = - let result = - Array.init (maxm+1) (fun _ -> - Array.init np (fun _ -> Array.make nq 0. ) ) - in - let empty = Array.make (maxm+1) 0. in - Array.iteri (fun ab shell_ab -> - let zero_m_array_tmp = - Array.mapi (fun cd shell_cd -> - if (abs_float coef.(ab).(cd) < cutoff) then - empty - else - let expo_pq_inv = - expo_inv_p.(ab) +. expo_inv_q.(cd) - in - let norm_pq_sq = - let x = (center_pq X).(ab).(cd) - and y = (center_pq Y).(ab).(cd) - and z = (center_pq Z).(ab).(cd) + let zm_array = Mat.init_cols np nq (fun i j -> + try + if (abs_float coef.{j,i} ) < 1.e-3*.cutoff then + raise NullQuartet; + + let expo_pq_inv = + expo_inv_p.{i} +. expo_inv_q.{j} in - x *. x +. y *. y +. z *. z - in - zero_m ~maxm ~expo_pq_inv ~norm_pq_sq - ) sq - in - (* Transpose result *) - let coef_ab = coef.(ab) in - for m=0 to maxm do - let result_m_ab = result.(m).(ab) - in - for cd=0 to nq-1 do - result_m_ab.(cd) <- zero_m_array_tmp.(cd).(m) *. coef_ab.(cd) - done - done - ) sp; - result - in + let center_pq = + Co.(Psp.center sp.(i-1) |- Psp.center sq.(j-1)) + in + let norm_pq_sq = + Co.dot center_pq center_pq + in - let norm = - let norm_coef_scale_q = - Csp.norm_scales shell_q + let zero_m_array = + zero_m ~maxm:0 ~expo_pq_inv ~norm_pq_sq + in + zero_m_array.(0) + with NullQuartet -> 0. + ) in + Mat.gemm_trace zm_array coef + with (Invalid_argument _) -> 0. + end + | _ -> + + let coef = + Array.init np (fun l -> Array.init nq (fun k -> cq.(k) *. cp.(l)) ) in - Array.to_list norm_coef_scale_p - |> List.map (fun v1 -> - Array.map (fun v2 -> v1 *. v2) norm_coef_scale_q) - |> Array.concat - in - let map_1d = Zmap.create (4*maxm) - and map_2d = Array.init (maxm+1) (fun _ -> Zmap.create (Array.length class_indices)) - in - (* Compute the integral class from the primitive shell quartet *) - Array.iteri (fun i key -> - let (angMom_a,angMom_b,angMom_c,angMom_d) = - match Zkey.to_powers key with - | Zkey.Twelve x -> x - | _ -> assert false + let norm = Cspc.norm_scales shell_pair_couple in + + let expo_inv_p = + Array.map (fun shell_ab -> Psp.exponent_inv shell_ab) sp + and expo_inv_q = + Array.map (fun shell_cd -> Psp.exponent_inv shell_cd) sq + in + + let expo_b = + Array.map (fun shell_ab -> Ps.exponent (Psp.shell_b shell_ab) ) sp + and expo_d = + Array.map (fun shell_cd -> Ps.exponent (Psp.shell_b shell_cd) ) sq + in + + let center_pq = + let result = + Array.init 3 (fun xyz -> + Array.init np (fun ab -> + let shell_ab = sp.(ab) in + Array.init nq (fun cd -> + let shell_cd = sq.(cd) + in + let cpq = + Co.(Psp.center shell_ab |- Psp.center shell_cd) + in + match xyz with + | 0 -> Co.get X cpq; + | 1 -> Co.get Y cpq; + | _ -> Co.get Z cpq; + ) + ) + ) + in function + | Co.X -> result.(0) + | Co.Y -> result.(1) + | Co.Z -> result.(2) + in + let center_pa = + let result = + Array.init 3 (fun xyz -> + Array.init np (fun ab -> + let shell_ab = sp.(ab) in + let cpa = + Co.(Psp.center shell_ab |- Cs.center shell_a) + in + match xyz with + | 0 -> Co.(get X cpa); + | 1 -> Co.(get Y cpa); + | _ -> Co.(get Z cpa); + ) + ) + in function + | Co.X -> result.(0) + | Co.Y -> result.(1) + | Co.Z -> result.(2) + in + let center_qc = + let result = + Array.init 3 (fun xyz -> + Array.init nq (fun cd -> + let shell_cd = sq.(cd) in + let cqc = + Co.(Psp.center shell_cd |- Cs.center shell_c) + in + match xyz with + | 0 -> Co.(get X cqc); + | 1 -> Co.(get Y cqc); + | _ -> Co.(get Z cqc); + ) + ) + in function + | Co.X -> result.(0) + | Co.Y -> result.(1) + | Co.Z -> result.(2) + in + let zero_m_array = + let result = + Array.init (maxm+1) (fun _ -> + Array.init np (fun _ -> Array.make nq 0. ) ) + in + let empty = Array.make (maxm+1) 0. in + Array.iteri (fun ab shell_ab -> + let zero_m_array_tmp = + Array.mapi (fun cd shell_cd -> + if (abs_float coef.(ab).(cd) < cutoff) then + empty + else + let expo_pq_inv = + expo_inv_p.(ab) +. expo_inv_q.(cd) + in + let norm_pq_sq = + let x = (center_pq X).(ab).(cd) + and y = (center_pq Y).(ab).(cd) + and z = (center_pq Z).(ab).(cd) + in + x *. x +. y *. y +. z *. z + in + + zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + ) sq in - try - if monocentric then - begin + (* Transpose result *) + let coef_ab = coef.(ab) in + for m=0 to maxm do + let result_m_ab = result.(m).(ab) + in + for cd=0 to nq-1 do + result_m_ab.(cd) <- zero_m_array_tmp.(cd).(m) *. coef_ab.(cd) + done + done + ) sp; + result + in + + let map_1d = Zmap.create (4*maxm) + and map_2d = Array.init (maxm+1) (fun _ -> Zmap.create (Array.length class_indices)) + in + (* Compute the integral class from the primitive shell quartet *) + Array.iteri (fun i key -> + let (angMom_a,angMom_b,angMom_c,angMom_d) = + match Zkey.to_powers key with + | Zkey.Twelve x -> x + | _ -> assert false + in + try + if monocentric then + begin if ( ((1 land angMom_a.x + angMom_b.x + angMom_c.x + angMom_d.x) =1) || ((1 land angMom_a.y + angMom_b.y + angMom_c.y + angMom_d.y) =1) || ((1 land angMom_a.z + angMom_b.z + angMom_c.z + angMom_d.z) =1) - ) then - raise NullQuartet - end; + ) then + raise NullQuartet + end; - (* Schwartz screening *) - if (np+nq> 24) then + (* Schwartz screening *) + if (np+nq> 24) then ( let schwartz_p = let key = Zkey.of_powers_twelve - angMom_a angMom_b angMom_a angMom_b + angMom_a angMom_b angMom_a angMom_b in match schwartz_p with | None -> 1. @@ -822,36 +789,37 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q if schwartz_p < cutoff then raise NullQuartet; let schwartz_q = let key = Zkey.of_powers_twelve - angMom_c angMom_d angMom_c angMom_d + angMom_c angMom_d angMom_c angMom_d in match schwartz_q with | None -> 1. | Some schwartz_q -> Zmap.find schwartz_q key in if schwartz_p *. schwartz_q < cutoff2 then raise NullQuartet; - ); - - let integral = - hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) - zero_m_array - (expo_b, expo_d) - (expo_inv_p, expo_inv_q) - (Csp.center_ab shell_p, - Csp.center_ab shell_q, center_pq) - (center_pa, center_qc) - map_1d map_2d np nq - in - contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i) - with NullQuartet -> () - ) class_indices + ); - end; + let integral = + hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) + zero_m_array + (expo_b, expo_d) + (expo_inv_p, expo_inv_q) + (Csp.a_minus_b shell_p, + Csp.a_minus_b shell_q, center_pq) + (center_pa, center_qc) + map_1d map_2d np nq + in + contracted_class.(i) <- contracted_class.(i) +. integral *. norm.(i) + with NullQuartet -> () + ) class_indices - let result = - Zmap.create (Array.length contracted_class) - in - Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; - result + end; + + let result = + Zmap.create (Array.length contracted_class) + in + Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; + result + with NullQuartet -> empty @@ -864,5 +832,5 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = match shell_p, shell_q with | Some shell_p, Some shell_q -> contracted_class_shell_pairs ~zero_m shell_p shell_q - | _ -> Zmap.create 0 + | _ -> empty