diff --git a/Basis/Basis.ml b/Basis/Basis.ml index a56f7c7..2ad62bd 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -28,13 +28,13 @@ let of_nuclei_and_general_basis n b = if (i > 0) then result.(i) <- Cs.with_index x ( result.(i-1).index + - (Array.length result.(i-1).powers )) + (Array.length result.(i-1).norm_coef_scale )) ) result ; let size = let n = ref 0 in for i=0 to (Array.length result) - 1 do - n := !n + (Array.length (result.(i).powers )) + n := !n + (Array.length (result.(i).norm_coef_scale)) done; !n in { contracted_shells = result ; size } diff --git a/Basis/ContractedShell.ml b/Basis/ContractedShell.ml index d7d08f5..168958c 100644 --- a/Basis/ContractedShell.ml +++ b/Basis/ContractedShell.ml @@ -10,19 +10,12 @@ type t = { size : int; norm_coef : float array; norm_coef_scale : float array; - powers : Zkey.t array; index : int; } module Am = AngularMomentum -(** Normalization coefficient of contracted function i, which depends on the - exponent and the angular momentum. Two conventions can be chosen : a single - normalisation factor for all functions of the class, or a coefficient which - depends on the powers of x,y and z. - Returns, for each contracted function, an array of functions taking as - argument the [|x;y;z|] powers. -*) + let compute_norm_coef expo totAngMom = let atot = Am.to_int totAngMom @@ -67,7 +60,7 @@ let make ~index ~expo ~coef ~center ~totAngMom = ) powers in { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ; - norm_coef_scale ; powers } + norm_coef_scale } let with_index a i = @@ -79,7 +72,7 @@ let to_string s = let open Printf in (match s.totAngMom with | Am.S -> sprintf "%3d " (s.index+1) - | _ -> sprintf "%3d-%-3d" (s.index+1) (s.index+(Array.length s.powers)) + | _ -> sprintf "%3d-%-3d" (s.index+1) (s.index+(Array.length s.norm_coef_scale)) ) ^ ( sprintf "%1s %8.3f %8.3f %8.3f " (Am.to_string s.totAngMom) (get X coord) (get Y coord) (get Z coord) ) ^ diff --git a/Basis/ContractedShell.mli b/Basis/ContractedShell.mli index c09cc4e..cacf36c 100644 --- a/Basis/ContractedShell.mli +++ b/Basis/ContractedShell.mli @@ -35,8 +35,7 @@ type t = private { totAngMom : AngularMomentum.t; (** Total angular momentum : {% $l = n_x + n_y + n_z$ %} *) size : int; (** Number of contracted functions, {% $m$ %} in the formula *) norm_coef : float array; (** Normalization coefficients of primitive functions {% $\mathcal{N}_i$ %} *) - norm_coef_scale : float array; (** Scaling factors {% $f_i$ %}, given in the same order as [powers]. *) - powers : Zkey.t array; (** Triplets {% $(n_x,n_y,n_z)$ %} encoded in a {!Zkey.t}. *) + norm_coef_scale : float array; (** Scaling factors {% $f_i$ %}, given in the same order as [AngularMomentum.zkey_array totAngMom]. *) index : int; (** Index in the basis set, represented as an array of contracted shells. *) } diff --git a/Basis/ERI.ml b/Basis/ERI.ml index ccae1b8..81b0ade 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -6,6 +6,7 @@ open Bigarray type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t +module Am = AngularMomentum module Bs = Basis module Cs = ContractedShell module Csp = ContractedShellPair @@ -212,10 +213,10 @@ let of_basis basis = ) else out := !out + 1; - ) shell.(l).powers - ) shell.(k).powers - ) shell.(j).powers - ) shell.(i).powers + ) Am.(zkey_array (Singlet shell.(l).Cs.totAngMom)) + ) Am.(zkey_array (Singlet shell.(k).Cs.totAngMom)) + ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) + ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) with NullIntegral -> () done; done; diff --git a/Basis/GamessReader.ml b/Basis/GamessReader.ml index a7e5c38..4b17765 100644 --- a/Basis/GamessReader.ml +++ b/Basis/GamessReader.ml @@ -1,3 +1,4 @@ +(** Reads a basis set in GAMESS format *) let read_basis filename = let lexbuf = let ic = open_in filename in diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index 3105162..a65be7b 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -155,8 +155,8 @@ let of_basis basis = in result.{i_c,j_c} <- value; result.{j_c,i_c} <- value; - ) shell.(i).powers - ) shell.(j).powers + ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) + ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) done; done; Mat.detri result; diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml index cf3d472..24ca506 100644 --- a/Basis/NucInt.ml +++ b/Basis/NucInt.ml @@ -6,7 +6,9 @@ open Lacaml.D type t = Mat.t +module Am = AngularMomentum module Bs = Basis +module Cs = ContractedShell (** (0|0)^m : Fundamental electron-nucleus repulsion integral $ \int \phi_p(r1) 1/r_{C} dr_1 $ @@ -86,8 +88,8 @@ let of_basis_nuclei basis nuclei = in eni_array.{j_c,i_c} <- value; eni_array.{i_c,j_c} <- value; - ) shell.(j).powers - ) shell.(i).powers + ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) + ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) done; done; Mat.detri eni_array; diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index bced594..48135ee 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -127,8 +127,8 @@ let of_basis basis = in result.{i_c,j_c} <- value; result.{j_c,i_c} <- value; - ) shell.(i).powers - ) shell.(j).powers + ) Am.(zkey_array (Singlet shell.(i).Cs.totAngMom)) + ) Am.(zkey_array (Singlet shell.(j).Cs.totAngMom)) done; done; Mat.detri result; diff --git a/Basis/Overlap_primitives.ml b/Basis/Overlap_primitives.ml index ea6dbfc..41eb0c3 100644 --- a/Basis/Overlap_primitives.ml +++ b/Basis/Overlap_primitives.ml @@ -8,22 +8,23 @@ let hvrr (angMom_a, angMom_b) (expo_inv_p) (center_ab, center_pa) (** Vertical recurrence relations *) let rec vrr angMom_a = - if angMom_a < 0 then 0. - else if angMom_a = 0 then 1. - else - (chop center_pa (fun () -> vrr (angMom_a-1))) - +. chop (0.5 *. (float_of_int (angMom_a-1)) *. expo_inv_p) + if angMom_a > 0 then + chop center_pa (fun () -> vrr (angMom_a-1)) + +. chop (0.5 *. float_of_int (angMom_a-1) *. expo_inv_p) (fun () -> vrr (angMom_a-2)) + else if angMom_a = 0 then 1. + else 0. (** Horizontal recurrence relations *) and hrr angMom_a angMom_b = - if angMom_b < 0 then 0. - else if angMom_b = 0 then vrr angMom_a - else + if angMom_b > 0 then hrr (angMom_a+1) (angMom_b-1) +. chop center_ab (fun () -> hrr angMom_a (angMom_b-1) ) + else if angMom_b = 0 then + vrr angMom_a + else 0. in hrr angMom_a angMom_b diff --git a/HartreeFock/Fock.ml b/HartreeFock/Fock.ml index 94b4481..b62184d 100644 --- a/HartreeFock/Fock.ml +++ b/HartreeFock/Fock.ml @@ -1,5 +1,6 @@ open Lacaml.D open Simulation +open Constants type t = Mat.t @@ -17,10 +18,11 @@ let make ~density simulation = for nu = 1 to nBas do for lambda = 1 to nBas do let p = m_P.{lambda,sigma} in - for mu = 1 to nu do - m_F.{mu,nu} <- m_F.{mu,nu} +. p *. - (m_G.{mu,lambda,nu,sigma} -. 0.5 *. m_G.{mu,lambda,sigma,nu}) - done + if abs_float p > epsilon then + for mu = 1 to nu do + m_F.{mu,nu} <- m_F.{mu,nu} +. p *. + (m_G.{mu,lambda,nu,sigma} -. 0.5 *. m_G.{mu,lambda,sigma,nu}) + done done done done; diff --git a/Utils/AngularMomentum.ml b/Utils/AngularMomentum.ml index bea478d..fea8c85 100644 --- a/Utils/AngularMomentum.ml +++ b/Utils/AngularMomentum.ml @@ -2,7 +2,7 @@ open Powers exception AngularMomentumError of string -type t = +type t = | S | P | D | F | G | H | I | J | K | L | M | N | O let of_char = function @@ -68,8 +68,12 @@ let n_functions a = (a*a + 3*a + 2)/2 +let zkey_array_memo : (kind, Zkey.t array) Hashtbl.t = + Hashtbl.create 13 + (** Returns an array of Zkeys corresponding to all possible angular momenta *) let zkey_array a = + let keys_1d l = let create_z { x ; y ; _ } = Powers.of_int_tuple (x,y,l-(x+y)) @@ -93,41 +97,50 @@ let zkey_array a = |> List.rev in - begin - match a with - | Singlet l1 -> - List.map (fun x -> Zkey.of_powers_three x) (keys_1d @@ to_int l1) + try + Hashtbl.find zkey_array_memo a - | Doublet (l1, l2) -> - List.map (fun a -> - List.map (fun b -> Zkey.of_powers_six a b) (keys_1d @@ to_int l2) - ) (keys_1d @@ to_int l1) - |> List.concat + with Not_found -> - | Triplet (l1, l2, l3) -> + let result = + begin + match a with + | Singlet l1 -> + List.map (fun x -> Zkey.of_powers_three x) (keys_1d @@ to_int l1) - List.map (fun a -> - List.map (fun b -> - List.map (fun c -> - Zkey.of_powers_nine a b c) (keys_1d @@ to_int l3) - ) (keys_1d @@ to_int l2) + | Doublet (l1, l2) -> + List.map (fun a -> + List.map (fun b -> Zkey.of_powers_six a b) (keys_1d @@ to_int l2) + ) (keys_1d @@ to_int l1) |> List.concat - ) (keys_1d @@ to_int l1) - |> List.concat - | Quartet (l1, l2, l3, l4) -> + | Triplet (l1, l2, l3) -> - List.map (fun a -> - List.map (fun b -> - List.map (fun c -> - List.map (fun d -> - Zkey.of_powers_twelve a b c d) (keys_1d @@ to_int l4) - ) (keys_1d @@ to_int l3) + List.map (fun a -> + List.map (fun b -> + List.map (fun c -> + Zkey.of_powers_nine a b c) (keys_1d @@ to_int l3) + ) (keys_1d @@ to_int l2) |> List.concat - ) (keys_1d @@ to_int l2) + ) (keys_1d @@ to_int l1) |> List.concat - ) (keys_1d @@ to_int l1) - |> List.concat - end - |> Array.of_list + + | Quartet (l1, l2, l3, l4) -> + + List.map (fun a -> + List.map (fun b -> + List.map (fun c -> + List.map (fun d -> + Zkey.of_powers_twelve a b c d) (keys_1d @@ to_int l4) + ) (keys_1d @@ to_int l3) + |> List.concat + ) (keys_1d @@ to_int l2) + |> List.concat + ) (keys_1d @@ to_int l1) + |> List.concat + end + |> Array.of_list + in + Hashtbl.add zkey_array_memo a result; + result diff --git a/Utils/NonNegativeFloat.ml b/Utils/NonNegativeFloat.ml index 2ee5ecc..43391c2 100644 --- a/Utils/NonNegativeFloat.ml +++ b/Utils/NonNegativeFloat.ml @@ -1,4 +1,4 @@ -type t = float +type t = float let of_float x = if x < 0. then invalid_arg (__FILE__^": of_float"); diff --git a/Utils/Zkey.ml b/Utils/Zkey.ml index 43acebd..9e5f728 100644 --- a/Utils/Zkey.ml +++ b/Utils/Zkey.ml @@ -2,8 +2,8 @@ open Powers type t = { - mutable left : int ; - mutable right : int ; + mutable left : int; + mutable right : int; kind : int ; } diff --git a/Utils/Zkey.mli b/Utils/Zkey.mli index 4c636b0..ec43a41 100644 --- a/Utils/Zkey.mli +++ b/Utils/Zkey.mli @@ -67,3 +67,4 @@ val equal : t -> t -> bool val compare : t -> t -> int (** Comparison function, used for sorting. *) +