diff --git a/Basis/Basis.ml b/Basis/Basis.ml index bc3730c..7af8434 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -37,3 +37,4 @@ Momentum X Y Z |> String.concat line) ^ line + diff --git a/Basis/ERI.ml b/Basis/ERI.ml index f88e5f0..1d20e0f 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -1,9 +1,6 @@ open Util -let cutoff = 1.e-20 -let log_cutoff = -. (log cutoff) - -(** (00|00)^m : Fundamental integral +(** (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 @@ -24,260 +21,7 @@ let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = ) -(** 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, totAngMom shell_b, - totAngMom shell_c, 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 - - +let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t = + TwoElectronRR.contracted_class ~zero_m shell_a shell_b shell_c shell_d diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml new file mode 100644 index 0000000..1ca3367 --- /dev/null +++ b/Basis/TwoElectronRR.ml @@ -0,0 +1,260 @@ +open Util + +(** 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 + + + + + +(** Computes all the two-electron integrals of the contracted shell quartet *) +let contracted_class ~zero_m 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, totAngMom shell_b, + totAngMom shell_c, 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/Utils/Util.ml b/Utils/Util.ml index 6fb0b6f..ae14763 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -1,3 +1,5 @@ +let cutoff = 1.e-20 + (** Constants *) let pi = acos (-1.) let pi_inv = 1. /. pi