From 4a4809fc7ec00d00a8e6b0190b605673eb0f61dc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 22 Jan 2018 23:19:24 +0100 Subject: [PATCH] Inserted overlap --- Basis/Basis.ml | 17 ++++++ Basis/Basis.mli | 4 ++ Basis/OneElectronRR.ml | 13 ++--- Basis/Overlap.ml | 111 ++++++++++++++++++++++++++++++++++++ Basis/Overlap_primitives.ml | 34 +++++++++++ Basis/Shell_pair.ml | 14 ++++- Basis/TwoElectronRR.ml | 5 -- Nuclei/Nuclei.ml | 11 ++++ Utils/Angular_momentum.ml | 4 +- Utils/Util.ml | 5 ++ Utils/Zkey.ml | 45 ++++++--------- run_integrals.ml | 46 ++++----------- 12 files changed, 230 insertions(+), 79 deletions(-) create mode 100644 Basis/Overlap.ml create mode 100644 Basis/Overlap_primitives.ml diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 5f9ddc1..7c1b869 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -45,4 +45,21 @@ let to_string b = ) ^ line +let file : string option ref = ref None + +let set_file f = + file := Some f + +let general_basis = lazy( + match !file with + | None -> failwith "basis set file not defined" + | Some filename -> Gamess_reader.read ~filename +) + + +let basis = lazy ( + of_nuclei_and_general_basis + (Lazy.force Nuclei.nuclei) (Lazy.force general_basis) +) + diff --git a/Basis/Basis.mli b/Basis/Basis.mli index ac72149..9626e93 100644 --- a/Basis/Basis.mli +++ b/Basis/Basis.mli @@ -1,3 +1,7 @@ type t = Contracted_shell.t array val of_nuclei_and_general_basis : Nuclei.t -> General_basis.t list -> t val to_string : t -> string +val set_file : string -> unit +val general_basis : + (Element.t * General_basis.general_contracted_shell array) list lazy_t +val basis : Contracted_shell.t array lazy_t diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index 2719b61..31aea67 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -129,27 +129,22 @@ let contracted_class_nuc ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t in () -(* TODO (* 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 coef_prod = - shell_p.(ab).Shell_pair.coef *. shell_q.(cd).Shell_pair.coef + shell_p.(ab).Shell_pair.coef in (** Screening on thr product of coefficients *) if (abs_float coef_prod) > 1.e-4*.cutoff then begin let expo_pq_inv = - shell_p.(ab).Shell_pair.expo_inv +. shell_q.(cd).Shell_pair.expo_inv + shell_p.(ab).Shell_pair.expo_inv in - let center_pq = - Coordinate.(shell_p.(ab).Shell_pair.center |- shell_q.(cd).Shell_pair.center) + let center_ab = + Coordinate.shell_p.(ab).Shell_pair.center_ab in let norm_pq_sq = Coordinate.dot center_pq center_pq diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml new file mode 100644 index 0000000..5c8b2cb --- /dev/null +++ b/Basis/Overlap.ml @@ -0,0 +1,111 @@ +open Util + +(** Computes all the overlap integrals of the contracted shell pair *) +let contracted_class shell_a shell_b : float Zmap.t = + + let shell_p = + Shell_pair.create_array shell_a shell_b + in + + (* Pre-computation of integral class indices *) + let class_indices = + Angular_momentum.zkey_array (Angular_momentum.Doublet + Contracted_shell.(totAngMom shell_a, totAngMom shell_b)) + 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 coef_prod = + shell_p.(ab).Shell_pair.coef + in + (** Screening on thr product of coefficients *) + if (abs_float coef_prod) > 1.e-4*.cutoff then + begin + let center_ab = + shell_p.(ab).Shell_pair.center_ab + in + let center_a = + shell_p.(ab).Shell_pair.center_a + in + let expo_inv = + shell_p.(ab).Shell_pair.expo_inv + in + + Array.iteri (fun i key -> + let (angMomA,angMomB) = + let a = Zkey.to_int_array Zkey.Kind_6 key in + ( [| a.(0) ; a.(1) ; a.(2) |], + [| a.(3) ; a.(4) ; a.(5) |] ) + in + let norm = + shell_p.(ab).Shell_pair.norm_fun angMomA angMomB + in + let f k = + Overlap_primitives.hvrr (angMomA.(k), angMomB.(k)) + expo_inv + (Coordinate.coord center_ab k, + Coordinate.coord center_a k) + in + let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in + contracted_class.(i) <- contracted_class.(i) +. coef_prod *. integral + ) class_indices + end + 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 + + + + + +(** Write all overlap integrals to a file *) +let to_file ~filename basis = + let to_int_tuple x = + let open Zkey in + match to_int_tuple Kind_3 x with + | Three x -> x + | _ -> assert false + in + + let oc = open_out filename in + for i=0 to (Array.length basis) - 1 do + print_int basis.(i).Contracted_shell.indice ; print_newline (); + for j=0 to i do + (* Compute all the integrals of the class *) + let cls = + contracted_class basis.(i) basis.(j) + in + + (* Write the data in the output file *) + Array.iteri (fun i_c powers_i -> + let i_c = basis.(i).Contracted_shell.indice + i_c + 1 in + let xi = to_int_tuple powers_i in + Array.iteri (fun j_c powers_j -> + let j_c = basis.(j).Contracted_shell.indice + j_c + 1 in + let xj = to_int_tuple powers_j in + let key = + Zkey.of_int_tuple (Zkey.Six (xi,xj)) + in + let value = + try + Zmap.find cls key + with Not_found -> failwith "Bug in overlap integrals" + in + if (abs_float value > cutoff) then + Printf.fprintf oc "%4d %4d %20.12e\n" + i_c j_c value + ) basis.(j).Contracted_shell.powers + ) basis.(i).Contracted_shell.powers; + done; + done; + close_out oc + diff --git a/Basis/Overlap_primitives.ml b/Basis/Overlap_primitives.ml new file mode 100644 index 0000000..1fe7a24 --- /dev/null +++ b/Basis/Overlap_primitives.ml @@ -0,0 +1,34 @@ +open Util + + +(** Horizontal and Vertical Recurrence Relations (HVRR) *) +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) + (fun () -> vrr (angMom_a-2)) + + + (** 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 + hrr (angMom_a+1) (angMom_b-1) +. chop center_ab + (fun () -> hrr angMom_a (angMom_b-1) ) + + in + hrr angMom_a angMom_b + + + +let overlap = hvrr + diff --git a/Basis/Shell_pair.ml b/Basis/Shell_pair.ml index 7c46bf8..a770455 100644 --- a/Basis/Shell_pair.ml +++ b/Basis/Shell_pair.ml @@ -4,6 +4,7 @@ type t = { expo : float; expo_inv : float; center_ab: Coordinate.t; + center_a : Coordinate.t; center : Coordinate.t; norm_sq : float; norm : float; @@ -38,6 +39,12 @@ let create_array ?cutoff p_a p_b = Array.init (Contracted_shell.size p_b) (fun j -> try + if Contracted_shell.center p_a = Contracted_shell.center p_b && + ((Angular_momentum.to_int @@ Contracted_shell.totAngMom p_a) + + (Angular_momentum.to_int @@ Contracted_shell.totAngMom p_b)) land 1 = 1 + then + raise Null_contribution; + let f2 = Contracted_shell.norm_coef p_b j in @@ -76,7 +83,10 @@ let create_array ?cutoff p_a p_b = 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 } + let center_a = + Coordinate.(center |- Contracted_shell.center p_a) + in + Some { i ; j ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_a ; center_ab ; norm_sq } with | Null_contribution -> None ) @@ -88,3 +98,5 @@ let create_array ?cutoff p_a p_b = |> List.map (function Some x -> x | None -> assert false) |> Array.of_list +open Util + diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index e3740f2..3929e68 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -1,10 +1,5 @@ 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 hvrr_two_e m (angMom_a, angMom_b, angMom_c, angMom_d) diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index 5069939..e9d5712 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -49,3 +49,14 @@ let to_string atoms = " +let file : string option ref = ref None + +let set_file f = + file := Some f + +let nuclei = lazy( + match !file with + | None -> failwith "coordinate file not defined" + | Some filename -> of_xyz_file ~filename +) + diff --git a/Utils/Angular_momentum.ml b/Utils/Angular_momentum.ml index 728fa64..299619a 100644 --- a/Utils/Angular_momentum.ml +++ b/Utils/Angular_momentum.ml @@ -99,8 +99,8 @@ let zkey_array a = | Doublet (l1, l2) -> List.map (fun a -> List.map (fun b -> - Zkey.of_int_tuple (Zkey.Six (a,b))) (keys_1d @@ to_int l1) - ) (keys_1d @@ to_int l2) + Zkey.of_int_tuple (Zkey.Six (a,b))) (keys_1d @@ to_int l2) + ) (keys_1d @@ to_int l1) |> List.concat | Triplet (l1, l2, l3) -> diff --git a/Utils/Util.ml b/Utils/Util.ml index 8ebca5b..4c9e76e 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -69,3 +69,8 @@ let rec pow a = function ;; +(** 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 ()) + diff --git a/Utils/Zkey.ml b/Utils/Zkey.ml index bd2b915..77fe828 100644 --- a/Utils/Zkey.ml +++ b/Utils/Zkey.ml @@ -10,18 +10,18 @@ type kind_array = | Kind_2 | Kind_1 +let (<|) x a = + Z.logor (Z.shift_left x 64) a + +let (<<) x a = + Int64.logor (Int64.shift_left x 10) (Int64.of_int a) + +let (<+) x a = + Int64.logor (Int64.shift_left x 16) (Int64.of_int a) + (** Build a Zkey from an array or 1, 2, 3, 4, 6, 9, or 12 integers *) let of_int_array ~kind a = - let (<|) x a = - Z.logor (Z.shift_left x 64) a - in - let (<<) x a = - Int64.logor (Int64.shift_left x 10) (Int64.of_int a) - in - let (<+) x a = - Int64.logor (Int64.shift_left x 16) (Int64.of_int a) - in match kind with | Kind_3 -> (Int64.of_int a.(0)) << a.(1) << a.(2) |> Z.of_int64 | Kind_6 -> (Int64.of_int a.(0)) << a.(1) << a.(2) << a.(3) << a.(4) << a.(5) @@ -56,15 +56,6 @@ type kind = | Nine of ((int*int*int)*(int*int*int)*(int*int*int)) | Twelve of ((int*int*int)*(int*int*int)*(int*int*int)*(int*int*int)) -let (<|) x a = - Z.logor (Z.shift_left x 64) a - -let (<<) x a = - Int64.logor (Int64.shift_left x 10) (Int64.of_int a) - -let (<+) x a = - Int64.logor (Int64.shift_left x 16) (Int64.of_int a) - let of_int_tuple a = match a with | One (a) -> Z.of_int a @@ -143,16 +134,8 @@ let to_int_array ~kind a = | Kind_1 -> [| Z.to_int a |] -let to_string ~kind a = - "< " ^ ( Z.to_string a ) ^ " | " ^ ( - to_int_array kind a - |> Array.map string_of_int - |> Array.to_list - |> String.concat ", " - ) ^ " >" - -(** Transform the Zkey into an int array *) +(** Transform the Zkey into an int tuple *) let to_int_tuple ~kind a = match kind with | Kind_3 -> Three ( Z.to_int @@ Z.extract a 20 10 , @@ -202,6 +185,14 @@ let to_int_tuple ~kind a = include Z +let to_string ~kind a = + "< " ^ ( Z.to_string a ) ^ " | " ^ ( + to_int_array kind a + |> Array.map string_of_int + |> Array.to_list + |> String.concat ", " + ) ^ " >" + (* let debug () = let k2 = of_int_array Kind_2 [| 1 ; 2 |] diff --git a/run_integrals.ml b/run_integrals.ml index 3637b15..77be74a 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -1,60 +1,36 @@ -let basis_file : string option ref = ref None -let coord_file : string option ref = ref None let out_file : string option ref = ref None let speclist = [ - ( "-b" , Arg.String (fun x -> basis_file := Some x) , + ( "-b" , Arg.String Basis.set_file , "File containing the atomic basis set") ; - ( "-c" , Arg.String (fun x -> coord_file := Some x) , + ( "-c" , Arg.String Nuclei.set_file , "File containing the nuclear coordinates") ; ( "-o" , Arg.String (fun x -> out_file := Some x) , "Output file") ; ] -let run ~coord ~basis ~out = - let coord_file = - match coord with - | None -> raise (Invalid_argument "Coordinate file should be specified with -c") - | Some x -> x - and basis_file = - match basis with - | None -> raise (Invalid_argument "Basis set file should be specified with -b") - | Some x -> x - and out_file = +let run ~out = + let out_file = match out with | None -> raise (Invalid_argument "Output file should be specified with -o") | Some x -> x in - let nuclei = - Nuclei.of_xyz_file ~filename:coord_file - in + let nuclei = Lazy.force Nuclei.nuclei in print_endline @@ Nuclei.to_string nuclei; - let basis = - let general_basis = - Gamess_reader.read ~filename:basis_file - in - Basis.of_nuclei_and_general_basis nuclei general_basis - in + let basis = Lazy.force Basis.basis in print_endline @@ Basis.to_string basis; - ERI.to_file ~filename:out_file basis - + ERI.to_file ~filename:(out_file^".eri") basis +(* + Overlap.to_file ~filename:(out_file^".overlap") basis + *) let () = let usage_msg = "Available options:" in Arg.parse speclist (fun _ -> ()) usage_msg; - run ~coord:!coord_file ~basis:!basis_file ~out:!out_file - (* - try - with - | Invalid_argument e -> - begin - print_string "Error: " ; print_endline e; print_newline (); - Arg.usage speclist usage_msg - end - *) + run ~out:!out_file