From 97e4a0966410f285930decb942e2aa662f0f54ec Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 17 Jan 2018 19:09:57 +0100 Subject: [PATCH] Cleaning --- Basis/Basis.ml | 1 - Basis/Contracted_shell.ml | 6 + Basis/ERI.ml | 283 +++++++++++++++++++++++++++++++++++++ Basis/Shell_pair.ml | 144 +++++++++---------- Utils/Angular_momentum.mli | 5 + 5 files changed, 359 insertions(+), 80 deletions(-) create mode 100644 Basis/ERI.ml diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 6e8e131..252b96b 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -1,5 +1,4 @@ (** General basis set read from a file *) - type primitive = { exponent: float ; coefficient: float diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml index 3138d5c..5bc414c 100644 --- a/Basis/Contracted_shell.ml +++ b/Basis/Contracted_shell.ml @@ -9,6 +9,12 @@ type t = { norm_coef : (int array -> float) array; } +let size a = a.size +let expo a i = a.expo.(i) +let coef a i = a.coef.(i) +let center a = a.center +let totAngMom a = a.totAngMom +let norm_coef a i = a.norm_coef.(i) (** Normalization coefficient of contracted function i, which depends on the diff --git a/Basis/ERI.ml b/Basis/ERI.ml new file mode 100644 index 0000000..d27432e --- /dev/null +++ b/Basis/ERI.ml @@ -0,0 +1,283 @@ +open Util + +let cutoff = 1.e-20 +let log_cutoff = -. (log cutoff) + +(** (00|00)^m : Fundamental 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 + boys_function ~maxm t + |> Array.mapi (fun m fm -> + two_over_sq_pi *. (if m mod 2 = 0 then fm else -.fm) *. + (pow exp_pq m) *. (sqrt exp_pq) + ) + + +(** In chop f g, evaluate g only if f is non zero, and return f *. (g ()) *) +let chop f g = + if (abs_float f) < cutoff then 0. + else f *. (g ()) + + +(** Horizontal and Vertical Recurrence Relations (HVRR) *) +let ghvrr m (angMom_a, angMom_b, angMom_c, angMom_d) + (totAngMom_a, totAngMom_b, totAngMom_c, totAngMom_d) + (maxm, zero_m_array) + (expo_b, expo_d) + (expo_inv_p, expo_inv_q) + (center_ab, center_cd, center_pq) + map + = + + let totAngMom_a = Angular_momentum.to_int totAngMom_a + and totAngMom_b = Angular_momentum.to_int totAngMom_b + and totAngMom_c = Angular_momentum.to_int totAngMom_c + and totAngMom_d = Angular_momentum.to_int totAngMom_d + in + + (** Vertical recurrence relations *) + let rec gvrr m angMom_a angMom_c totAngMom_a totAngMom_c = + + if angMom_a.(0) < 0 || angMom_a.(1) < 0 || angMom_a.(2) < 0 + || angMom_c.(0) < 0 || angMom_c.(1) < 0 || angMom_c.(2) < 0 then 0. + else + match (totAngMom_a, totAngMom_c) with + | (0,0) -> zero_m_array.(m) + | (_,0) -> + + let key = [| angMom_a.(0)+1; angMom_a.(1)+1; angMom_a.(2)+1; |] + |> Zkey.(of_int_array ~kind:Kind_3) + in + + let (found, result) = + try (true, Zmap.find map.(m) key) with + | Not_found -> (false, + let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] + and amm = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] + and xyz = + match angMom_a with + | [|0;0;_|] -> 2 + | [|0;_;_|] -> 1 + | _ -> 0 + in + am.(xyz) <- am.(xyz) - 1; + amm.(xyz) <- amm.(xyz) - 2; + chop (-. expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz)) + (fun () -> gvrr m am angMom_c (totAngMom_a-1) totAngMom_c ) + +. chop (expo_inv_p *. (Coordinate.coord center_pq xyz)) + (fun () -> gvrr (m+1) am angMom_c (totAngMom_a-1) totAngMom_c ) + +. chop ((float_of_int am.(xyz)) *. expo_inv_p *. 0.5) + (fun () -> gvrr m amm angMom_c (totAngMom_a-2) totAngMom_c + +. chop expo_inv_p (fun () -> + gvrr (m+1) amm angMom_c (totAngMom_a-2) totAngMom_c) ) ) + in + if not found then + Zmap.add map.(m) key result; + result + + | (_,_) -> + + let key = [| angMom_a.(0)+1; angMom_a.(1)+1; angMom_a.(2)+1; + angMom_c.(0)+1; angMom_c.(1)+1; angMom_c.(2)+1; |] + |> Zkey.(of_int_array ~kind:Kind_6) + in + + let (found, result) = + try (true, Zmap.find map.(m) key) with + | Not_found -> (false, + let am = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] + and cm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] + and cmm = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] + and xyz = + match angMom_c with + | [|0;0;_|] -> 2 + | [|0;_;_|] -> 1 + | _ -> 0 + in + am.(xyz) <- am.(xyz) - 1; + cm.(xyz) <- cm.(xyz) - 1; + cmm.(xyz) <- cmm.(xyz) - 2; + chop (-. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz) ) + (fun () -> gvrr m angMom_a cm totAngMom_a (totAngMom_c-1) ) + -. chop (expo_inv_q *. (Coordinate.coord center_pq xyz)) + (fun () -> gvrr (m+1) angMom_a cm totAngMom_a (totAngMom_c-1) ) + +. chop ((float_of_int cm.(xyz)) *. expo_inv_q *. 0.5 ) + (fun () -> gvrr m angMom_a cmm totAngMom_a (totAngMom_c-2) + +. chop expo_inv_q + (fun () -> gvrr (m+1) angMom_a cmm totAngMom_a (totAngMom_c-2) ) ) + -. chop ((float_of_int angMom_a.(xyz)) *. expo_inv_p *. expo_inv_q *. 0.5 ) + (fun () -> gvrr (m+1) am cm (totAngMom_a-1) (totAngMom_c-1) ) ) + in + if not found then + Zmap.add map.(m) key result; + result + + + + + (** Horizontal recurrence relations *) + and ghrr m angMom_a angMom_b angMom_c angMom_d + totAngMom_a totAngMom_b totAngMom_c totAngMom_d = + + if angMom_b.(0) < 0 || angMom_b.(1) < 0 || angMom_b.(2) < 0 + || angMom_d.(0) < 0 || angMom_d.(1) < 0 || angMom_d.(2) < 0 then 0. + else + match (totAngMom_b, totAngMom_d) with + | (0,0) -> gvrr m angMom_a angMom_c totAngMom_a totAngMom_c + | (_,_) -> + + let key = [| angMom_a.(0)+1; angMom_a.(1)+1; angMom_a.(2)+1; + angMom_b.(0)+1; angMom_b.(1)+1; angMom_b.(2)+1; + angMom_c.(0)+1; angMom_c.(1)+1; angMom_c.(2)+1; + angMom_d.(0)+1; angMom_d.(1)+1; angMom_d.(2)+1; |] + |> Zkey.(of_int_array ~kind:Kind_12) + in + + let (found, result) = + try (true, Zmap.find map.(m) key) with + | Not_found -> (false, + begin + match totAngMom_d with + | 0 -> + let ap = [| angMom_a.(0) ; angMom_a.(1) ; angMom_a.(2) |] + and bm = [| angMom_b.(0) ; angMom_b.(1) ; angMom_b.(2) |] + and xyz = + match angMom_b with + | [|0;0;_|] -> 2 + | [|0;_;_|] -> 1 + | _ -> 0 + in + ap.(xyz) <- ap.(xyz) + 1; + bm.(xyz) <- bm.(xyz) - 1; + ghrr m ap bm angMom_c angMom_d (totAngMom_a+1) (totAngMom_b-1) + totAngMom_c totAngMom_d + +. chop (Coordinate.coord center_ab xyz) (fun () -> + ghrr m angMom_a bm angMom_c angMom_d totAngMom_a (totAngMom_b-1) + totAngMom_c totAngMom_d ) + | _ -> + let cp = [| angMom_c.(0) ; angMom_c.(1) ; angMom_c.(2) |] + and dm = [| angMom_d.(0) ; angMom_d.(1) ; angMom_d.(2) |] + and xyz = + match angMom_d with + | [|0;0;_|] -> 2 + | [|0;_;_|] -> 1 + | _ -> 0 + in + cp.(xyz) <- cp.(xyz) + 1; + dm.(xyz) <- dm.(xyz) - 1; + ghrr m angMom_a angMom_b cp dm totAngMom_a totAngMom_b + (totAngMom_c+1) (totAngMom_d-1) + +. chop (Coordinate.coord center_cd xyz) (fun () -> + ghrr m angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b + totAngMom_c (totAngMom_d-1) ) + end) + in + if not found then + Zmap.add map.(m) key result; + result + + in + ghrr m angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b + totAngMom_c totAngMom_d + + + + + +(** Electron-electron repulsion integral *) +let erint_contracted_class shell_a shell_b shell_c shell_d : float Zmap.t = + + let shell_p = Shell_pair.create_array shell_a shell_b + and shell_q = Shell_pair.create_array shell_c shell_d + and maxm = + let open Angular_momentum in + (to_int @@ Contracted_shell.totAngMom shell_a) + (to_int @@ Contracted_shell.totAngMom shell_b) + + (to_int @@ Contracted_shell.totAngMom shell_c) + (to_int @@ Contracted_shell.totAngMom shell_d) + in + + (* Pre-computation of integral class indices *) + let class_indices = + Angular_momentum.zkey_array + (Angular_momentum.Kind_4 + (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, + Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d)) + in + + let contracted_class = + Array.make (Array.length class_indices) 0.; + in + + (* Compute all integrals in the shell for each pair of significant shell pairs *) + + for ab=0 to (Array.length shell_p - 1) + do + let b = shell_p.(ab).Shell_pair.j in + + for cd=0 to (Array.length shell_q - 1) + do + + let d = shell_q.(cd).Shell_pair.j in + + let expo_pq_inv = + shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv + in + let center_pq = + Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center) + in + let norm_pq_sq = + Coordinate.dot center_pq center_pq + in + + let zero_m_array = + zero_m ~maxm ~expo_pq_inv ~norm_pq_sq + in + let map = Array.init maxm (fun _ -> Zmap.create 129) in + (* Compute the integral class from the primitive shell quartet *) + Array.iteri (fun i key -> + let (angMomA,angMomB,angMomC,angMomD) = + let a = Zkey.to_int_array Zkey.Kind_12 key in + ( [| a.(0) ; a.(1) ; a.(2) |], + [| a.(3) ; a.(4) ; a.(5) |], + [| a.(6) ; a.(7) ; a.(8) |], + [| a.(9) ; a.(10) ; a.(11) |] ) + in + let integral = + ghvrr 0 (angMomA, angMomB, angMomC, angMomD) + (Contracted_shell.totAngMom shell_a, Contracted_shell.totAngMom shell_b, + Contracted_shell.totAngMom shell_c, Contracted_shell.totAngMom shell_d) + (maxm, zero_m_array) + (Contracted_shell.expo shell_b b, Contracted_shell.expo shell_d d) + (shell_p.(ab).Shell_pair.expo_inv, shell_q.(cd).Shell_pair.expo_inv) + (shell_p.(ab).Shell_pair.center_ab, shell_q.(cd).Shell_pair.center_ab, center_pq) + map + in + let norm = + shell_p.(ab).Shell_pair.norm_fun angMomA angMomB *. shell_q.(cd).Shell_pair.norm_fun angMomC angMomD + in + let coef_prod = + shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef *. norm + in + contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral + ) class_indices + + done + done; + let result = + Zmap.create (Array.length contracted_class) + in + Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; + result + + + diff --git a/Basis/Shell_pair.ml b/Basis/Shell_pair.ml index d2fb9e5..87d97c6 100644 --- a/Basis/Shell_pair.ml +++ b/Basis/Shell_pair.ml @@ -3,8 +3,8 @@ open Util type t = { expo : float; expo_inv : float; - center_ab: float array; - center : float array; + center_ab: Coordinate.t; + center : Coordinate.t; norm_sq : float; norm : float; coef : float; @@ -21,83 +21,69 @@ let create_array ?(cutoff=0.) p_a p_b = else -. (log cutoff) in - let x_a = Coordinate.x p_a.Contracted_shell.center - and y_a = Coordinate.y p_a.Contracted_shell.center - and z_a = Coordinate.z p_a.Contracted_shell.center - and x_b = Coordinate.x p_b.Contracted_shell.center - and y_b = Coordinate.y p_b.Contracted_shell.center - and z_b = Coordinate.z p_b.Contracted_shell.center + let center_ab = + Coordinate.(Contracted_shell.center p_a |- Contracted_shell.center p_b) in - (* - match p_a.Contracted_shell.center, p_b.Contracted_shell.center with - | [|x_a; y_a; z_a|], [|x_b; y_b; z_b|] -> - *) - let center_ab = - Coordinate.(p_a.Contracted_shell.center |- p_b.Contracted_shell.center) - in - let norm_sq = - Coordinate.dot center_ab center_ab - in - Array.init p_a.Contracted_shell.size (fun i -> - let p_a_expo_center = - [| p_a.Contracted_shell.expo.(i) *. x_a ; p_a.Contracted_shell.expo.(i) *. y_a ; p_a.Contracted_shell.expo.(i) *. z_a |] - in + let norm_sq = + Coordinate.dot center_ab center_ab + in + Array.init (Contracted_shell.size p_a) (fun i -> + let p_a_expo_center = + Coordinate.(Contracted_shell.expo p_a i |. Contracted_shell.center p_a) + in + let f1 = + Contracted_shell.norm_coef p_a i + in - Array.init p_b.Contracted_shell.size (fun j -> - try - let f1 = - p_a.Contracted_shell.norm_coef.(i) - in - let f2 = - p_b.Contracted_shell.norm_coef.(j) - in - let norm_fun a b = - f1 a *. f2 b - in - let norm = - norm_fun - [| Angular_momentum.to_int p_a.Contracted_shell.totAngMom ; 0 ; 0 |] - [| Angular_momentum.to_int p_b.Contracted_shell.totAngMom ; 0 ; 0 |] - in - if (norm < cutoff) then - raise Null_contribution; - let p_b_expo_center = - [| p_b.Contracted_shell.expo.(j) *. x_b ; p_b.Contracted_shell.expo.(j) *. y_b ; p_b.Contracted_shell.expo.(j) *. z_b |] - in - let expo = p_a.Contracted_shell.expo.(i) +. p_b.Contracted_shell.expo.(j) in - let expo_inv = 1. /. expo in - let center = - [| (p_a_expo_center.(0) +. p_b_expo_center.(0)) *. expo_inv; - (p_a_expo_center.(1) +. p_b_expo_center.(1)) *. expo_inv; - (p_a_expo_center.(2) +. p_b_expo_center.(2)) *. expo_inv |] - in - let argexpo = - p_a.Contracted_shell.expo.(i) *. p_b.Contracted_shell.expo.(j) - *. norm_sq *. expo_inv - in - if (argexpo > log_cutoff) then - raise Null_contribution; - let g = - (pi *. expo_inv)**(1.5) *. exp(-. argexpo) - in - let norm_inv = 1./.norm in - let norm_fun a b = - norm_inv *. norm_fun a b - in - let coef = - norm *. p_a.Contracted_shell.coef.(i) *. p_b.Contracted_shell.coef.(j) *. g - in - if (abs_float coef < cutoff) then - raise Null_contribution; - Some { i ; j ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_ab=(Coordinate.to_float_array center_ab) ; norm_sq } - with - | Null_contribution -> None - ) - ) - |> Array.to_list - |> Array.concat - |> Array.to_list - |> List.filter (function Some _ -> true | None -> false) - |> List.map (function Some x -> x | None -> assert false) - |> Array.of_list + Array.init (Contracted_shell.size p_b) (fun j -> + try + let f2 = + Contracted_shell.norm_coef p_b j + in + let norm_fun a b = + f1 a *. f2 b + in + let norm = + norm_fun + [| Angular_momentum.to_int @@ Contracted_shell.totAngMom p_a ; 0 ; 0 |] + [| Angular_momentum.to_int @@ Contracted_shell.totAngMom p_b ; 0 ; 0 |] + in + if (norm < cutoff) then + raise Null_contribution; + let p_b_expo_center = + Coordinate.(Contracted_shell.expo p_b j |. Contracted_shell.center p_b) + in + let expo = Contracted_shell.(expo p_a i +. expo p_b j) in + let expo_inv = 1. /. expo in + let center = + Coordinate.( expo_inv |. (p_a_expo_center |+ p_b_expo_center ) ) + in + let argexpo = + Contracted_shell.(expo p_a i *. expo p_b j) *. norm_sq *. expo_inv + in + if (argexpo > log_cutoff) then + raise Null_contribution; + let g = + (pi *. expo_inv)**(1.5) *. exp(-. argexpo) + in + let norm_inv = 1./.norm in + let norm_fun a b = + norm_inv *. norm_fun a b + in + let coef = + norm *. Contracted_shell.(coef p_a i *. coef p_b j) *. g + in + if (abs_float coef < cutoff) then + raise Null_contribution; + Some { i ; j ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_ab ; norm_sq } + with + | Null_contribution -> None + ) + ) + |> Array.to_list + |> Array.concat + |> Array.to_list + |> List.filter (function Some _ -> true | None -> false) + |> List.map (function Some x -> x | None -> assert false) + |> Array.of_list diff --git a/Utils/Angular_momentum.mli b/Utils/Angular_momentum.mli index c13b6bc..f410c38 100644 --- a/Utils/Angular_momentum.mli +++ b/Utils/Angular_momentum.mli @@ -2,3 +2,8 @@ exception AngularMomentumError of string type t = S | P | D | F | G | H | I | J | K | L | M | N | O val of_char : char -> t val to_string : t -> string +val to_char : t -> char +val to_int : t -> int +val of_int : int -> t +type kind = Kind_2 of (t * t) | Kind_4 of (t * t * t * t) +val zkey_array : kind -> Z.t array