From 73c8e48731df3931815d4d9c2fbfc007471dce40 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 19 Jan 2018 03:14:06 +0100 Subject: [PATCH] Optimized ERIs --- Basis/Basis.ml | 14 +-- Basis/Basis.mli | 3 + Basis/ERI.ml | 188 +++++++++++++++++++++++++++++++++++-- Basis/Gamess_parser.mly | 2 +- Basis/Gamess_reader.ml | 5 +- Basis/General_basis.ml | 4 +- Basis/TwoElectronRR.ml | 1 - Makefile | 1 + Nuclei/Nuclei.ml | 2 +- Utils/Angular_momentum.ml | 19 +++- Utils/Angular_momentum.mli | 3 +- Utils/Util.ml | 2 +- run_integrals.ml | 42 ++++++--- 13 files changed, 246 insertions(+), 40 deletions(-) create mode 100644 Basis/Basis.mli diff --git a/Basis/Basis.ml b/Basis/Basis.ml index 7af8434..bef0162 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -1,4 +1,4 @@ -type t +type t = Contracted_shell.t array (** Returns an array of the basis set per atom *) let of_nuclei_and_general_basis n b = @@ -12,6 +12,8 @@ let of_nuclei_and_general_basis n b = in Contracted_shell.create ~expo ~coef ~totAngMom ~center) ) n + |> Array.to_list + |> Array.concat let to_string b = @@ -27,14 +29,12 @@ Angular Coordinates (Bohr) Exponents Coefficients Momentum X Y Z ----------------------------------------------------------------------- " - ^( Array.mapi (fun atom_id basis -> - Array.map (fun i -> - Contracted_shell.to_string i) basis + ^ + ( Array.map (fun i -> + Contracted_shell.to_string i) b |> Array.to_list |> String.concat line - ) b - |> Array.to_list - |> String.concat line) + ) ^ line diff --git a/Basis/Basis.mli b/Basis/Basis.mli new file mode 100644 index 0000000..ac72149 --- /dev/null +++ b/Basis/Basis.mli @@ -0,0 +1,3 @@ +type t = Contracted_shell.t array +val of_nuclei_and_general_basis : Nuclei.t -> General_basis.t list -> t +val to_string : t -> string diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 1d20e0f..467422c 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -1,11 +1,15 @@ +(** Electron-electron repulsion integrals *) + open Util (** (00|00)^m : Fundamental electron repulsion integral - $ \int \int \phi_p(r1) 1/r_{12} \phi_q(r2) dr_1 dr_2 $ + $ \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$ + 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 = @@ -21,7 +25,179 @@ let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = ) -(** Electron-electron repulsion integral *) -let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t = +(** Compute all the integrals of a contracted class *) +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 + +type n_cls = { n : int ; cls : Z.t array } +(** Write all integrals to a file *) +let to_file ~filename basis = + let oc = open_out filename in + let key_array = Array.make 12 0 in + let zkey = Array.map (fun b -> + let result = + Angular_momentum.(zkey_array (Kind_1 b.Contracted_shell.totAngMom)) + in + { n=Array.length result ; cls=result } + ) basis + in + let i_shift = ref 1 in + let j_shift = ref 1 in + let k_shift = ref 1 in + let l_shift = ref 1 in + for i=0 to (Array.length basis) - 1 do + print_int !i_shift ; print_newline (); + j_shift := 1; + for j=0 to i do + k_shift := 1; + for k=0 to i do + l_shift := 1; + for l=0 to k do + let cls = + contracted_class basis.(i) basis.(j) basis.(k) basis.(l) + in + + for i_c = 0 to zkey.(i).n - 1 do + let x = Zkey.(to_int_array Kind_3 zkey.(i).cls.(i_c)) in + key_array.(0) <- x.(0); + key_array.(1) <- x.(1); + key_array.(2) <- x.(2); + for j_c = 0 to zkey.(j).n - 1 do + let x = Zkey.(to_int_array Kind_3 zkey.(j).cls.(j_c)) in + key_array.(3) <- x.(0); + key_array.(4) <- x.(1); + key_array.(5) <- x.(2); + for k_c = 0 to zkey.(k).n - 1 do + let x = Zkey.(to_int_array Kind_3 zkey.(k).cls.(k_c)) in + key_array.(6) <- x.(0); + key_array.(7) <- x.(1); + key_array.(8) <- x.(2); + for l_c = 0 to zkey.(l).n - 1 do + let x = Zkey.(to_int_array Kind_3 zkey.(l).cls.(l_c)) in + key_array.( 9) <- x.(0); + key_array.(10) <- x.(1); + key_array.(11) <- x.(2); + let key = + Zkey.(of_int_array Kind_12 key_array) + in + let value = + Zmap.find cls key + in + if (abs_float value > cutoff) then + Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" + (!i_shift+i_c) (!j_shift+j_c) (!k_shift+k_c) (!l_shift+l_c) + value + done + done + done + done; + l_shift := !l_shift + zkey.(l).n + done; + k_shift := !k_shift + zkey.(k).n + done; + j_shift := !j_shift + zkey.(j).n + done; + i_shift := !i_shift + zkey.(i).n + done + ; + close_out oc + + +(* +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.Contracted_shell.totAngMom)) + in + { n=Array.length result ; cls=result } + ) basis + in + + let key_array = Array.make 12 0 in + let result = ref [] in + + let i_shift = ref 1 in + for i=0 to (Array.length basis) - 1 do + print_int !i_shift ; print_newline (); + let j_shift = ref 1 in + for j=0 to i do + let k_shift = ref 1 in + for k=0 to i do + let l_shift = ref 1 in + for l=0 to k do + let cls = + contracted_class basis.(i) basis.(j) basis.(k) basis.(l) + in + + for i_c = 0 to zkey.(i).n - 1 do + let x = Zkey.(to_int_array Kind_3 zkey.(i).cls.(i_c)) in + key_array.(0) <- x.(0); + key_array.(1) <- x.(1); + key_array.(2) <- x.(2); + for j_c = 0 to zkey.(j).n - 1 do + let x = Zkey.(to_int_array Kind_3 zkey.(j).cls.(j_c)) in + key_array.(3) <- x.(0); + key_array.(4) <- x.(1); + key_array.(5) <- x.(2); + for k_c = 0 to zkey.(k).n - 1 do + let x = Zkey.(to_int_array Kind_3 zkey.(k).cls.(k_c)) in + key_array.(6) <- x.(0); + key_array.(7) <- x.(1); + key_array.(8) <- x.(2); + for l_c = 0 to zkey.(l).n - 1 do + let x = Zkey.(to_int_array Kind_3 zkey.(l).cls.(l_c)) in + key_array.( 9) <- x.(0); + key_array.(10) <- x.(1); + key_array.(11) <- x.(2); + let key = + Zkey.(of_int_array Kind_12 key_array) + in + let value = + Zmap.find cls key + in + if (abs_float value > cutoff) then + let key = + Zkey.of_int_array Zkey.Kind_4 [| + (!i_shift+i_c);(!j_shift+j_c);(!k_shift+k_c);(!l_shift+l_c) + |] + in + result := (key, value) :: !result + done + done + done + done; + l_shift := !l_shift + zkey.(l).n + done; + k_shift := !k_shift + zkey.(k).n + done; + j_shift := !j_shift + zkey.(j).n + done; + i_shift := !i_shift + zkey.(i).n + done + ; + + let result = Array.of_list !result in + let result = Zmap.create (Array.length result) in + + for i=0 to !i_shift - 2 do + for j=0 to !i_shift - 2 do + for k=0 to !i_shift - 2 do + for l=0 to !i_shift - 2 do + let key = + Zkey.of_int_array Zkey.Kind_4 [| i;j;k;l |] + in + try + let value = + Zmap.find result key + in + Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" + i j k l value + with Not_found -> () + done + done + done + done; + close_out oc +*) diff --git a/Basis/Gamess_parser.mly b/Basis/Gamess_parser.mly index e15978a..082e4f8 100644 --- a/Basis/Gamess_parser.mly +++ b/Basis/Gamess_parser.mly @@ -26,7 +26,7 @@ basis: | element shell_array EOF { ($1, $2) } element: - | ELEMENT { $1 } + | ELEMENT { Element.of_string $1 } shell_array: | shell_list { Array.of_list @@ List.rev $1 } diff --git a/Basis/Gamess_reader.ml b/Basis/Gamess_reader.ml index f6f69c5..e2a2b7e 100644 --- a/Basis/Gamess_reader.ml +++ b/Basis/Gamess_reader.ml @@ -7,12 +7,9 @@ let read ~filename = in let rec aux accu = try - let element, basis = + let key, basis = Gamess_parser.input Basis_lexer.read_all lexbuf in - let key = - Element.of_string element - in aux ((key, basis)::accu) with | Parsing.Parse_error -> List.rev accu diff --git a/Basis/General_basis.ml b/Basis/General_basis.ml index 7ea7642..edb20d3 100644 --- a/Basis/General_basis.ml +++ b/Basis/General_basis.ml @@ -4,11 +4,9 @@ type primitive = { coefficient: float } -type element_name = string - type general_contracted_shell = Angular_momentum.t * (primitive array) -type t = element_name * (general_contracted_shell array) +type t = Element.t * (general_contracted_shell array) let string_of_primitive ?id prim = diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index 1ca3367..1d9289b 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -257,4 +257,3 @@ let contracted_class ~zero_m shell_a shell_b shell_c shell_d : float Zmap.t = - diff --git a/Makefile b/Makefile index 7b5d1c0..dee7889 100644 --- a/Makefile +++ b/Makefile @@ -4,6 +4,7 @@ INCLUDE_DIRS=Nuclei,Utils,Basis LIBS= PKGS= OCAMLCFLAGS="-g" +OCAMLCFLAGS="-unsafe -noassert -safe-string" OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) -Is $(INCLUDE_DIRS) MLLFILES=$(wildcard */*.mll) $(wildcard *.mll) MLYFILES=$(wildcard */*.mly) $(wildcard *.mly) diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index a683900..5069939 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -1,4 +1,4 @@ -type coordinates = ( (Element.t * Coordinate.t) array) +type t = (Element.t * Coordinate.t) array let of_xyz_file ~filename = let lexbuf = diff --git a/Utils/Angular_momentum.ml b/Utils/Angular_momentum.ml index bd48173..7b443d6 100644 --- a/Utils/Angular_momentum.ml +++ b/Utils/Angular_momentum.ml @@ -53,10 +53,18 @@ let of_int = function type kind = +| Kind_1 of t | Kind_2 of (t*t) | Kind_4 of (t*t*t*t) +let n_functions a = + let a = + to_int a + in + (a*a + 3*a + 2)/2 + + (** Returns an array of Zkeys corresponding to all possible angular momenta *) let zkey_array a = let keys_1d l = @@ -84,6 +92,14 @@ let zkey_array a = begin match a with + | Kind_1 l1 -> + + let a = Array.init 3 (fun _ -> 0) in + List.map (fun (cx,cy,cz) -> + a.(0) <- cx ; a.(1) <- cy ; a.(2) <- cz ; + Zkey.(of_int_array Kind_3 a) + ) (keys_1d @@ to_int l1) + | Kind_2 (l1, l2) -> let a = Array.init 6 (fun _ -> 0) in @@ -94,6 +110,7 @@ let zkey_array a = Zkey.(of_int_array Kind_6 a) ) (keys_1d @@ to_int l1) ) (keys_1d @@ to_int l2) + |> List.concat | Kind_4 (l1, l2, l3, l4) -> @@ -113,7 +130,7 @@ let zkey_array a = ) (keys_1d @@ to_int l2) |> List.concat ) (keys_1d @@ to_int l1) + |> List.concat end - |> List.concat |> Array.of_list diff --git a/Utils/Angular_momentum.mli b/Utils/Angular_momentum.mli index f410c38..73576a6 100644 --- a/Utils/Angular_momentum.mli +++ b/Utils/Angular_momentum.mli @@ -5,5 +5,6 @@ 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) +type kind = Kind_1 of t | Kind_2 of (t * t) | Kind_4 of (t * t * t * t) +val n_functions : t -> int val zkey_array : kind -> Z.t array diff --git a/Utils/Util.ml b/Utils/Util.ml index ae14763..5a666ff 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -1,4 +1,4 @@ -let cutoff = 1.e-20 +let cutoff = 1.e-50 (** Constants *) let pi = acos (-1.) diff --git a/run_integrals.ml b/run_integrals.ml index 50f9d8d..2d04c56 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -1,44 +1,58 @@ 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) , - "File containing the atomic basis set.") ; -( "-c" , Arg.String (fun x -> coord_file := Some x) , - "File containing the nuclear coordinates.") ; + ( "-b" , Arg.String (fun x -> basis_file := Some x) , + "File containing the atomic basis set") ; + ( "-c" , Arg.String (fun x -> coord_file := Some x) , + "File containing the nuclear coordinates") ; + ( "-o" , Arg.String (fun x -> out_file := Some x) , + "Output file") ; ] -let run ~coord ~basis = +let run ~coord ~basis ~out = let coord_file = match coord with - | None -> raise (Invalid_argument "Coordinate file should be specified") + | 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") + | None -> raise (Invalid_argument "Basis set file should be specified with -b") + | Some x -> x + and 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 - and general_basis = - Gamess_reader.read ~filename:basis_file 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 - Basis.to_string basis - |> print_endline + print_endline @@ Basis.to_string basis; + + ERI.to_file ~filename:out_file basis + let () = let usage_msg = "Available options:" in Arg.parse speclist (fun _ -> ()) usage_msg; try - run ~coord:!coord_file ~basis:!basis_file + run ~coord:!coord_file ~basis:!basis_file ~out:!out_file with | Invalid_argument e -> - (print_string "Error: " ; print_endline e; print_newline () ; Arg.usage speclist usage_msg) - ; + begin + print_string "Error: " ; print_endline e; print_newline (); + Arg.usage speclist usage_msg + end +