From 11151c997e3d437163723d217601e2c7b5f676d1 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 16 Apr 2020 19:49:23 +0200 Subject: [PATCH] Introduced multipole integrals --- Basis/AOBasis.ml | 8 +- Basis/AOBasis.mli | 3 + Basis/Basis.ml | 6 ++ Basis/Basis.mli | 5 ++ Basis/GeneralBasis.ml | 40 ++++++--- Basis/GeneralBasis.mli | 16 ++-- Basis/Multipole.ml | 199 +++++++++++++++++++++++++++++++++++++++++ Basis/Multipole.mli | 34 +++++++ Basis/Overlap.ml | 124 ++++++++++++------------- Nuclei/Nuclei.ml | 30 +++++-- Nuclei/Nuclei.mli | 3 + 11 files changed, 383 insertions(+), 85 deletions(-) create mode 100644 Basis/Multipole.ml create mode 100644 Basis/Multipole.mli diff --git a/Basis/AOBasis.ml b/Basis/AOBasis.ml index fa2cd0b..4499393 100644 --- a/Basis/AOBasis.ml +++ b/Basis/AOBasis.ml @@ -5,6 +5,7 @@ type t = { basis : Basis.t ; overlap : Overlap.t lazy_t; + multipole : Multipole.t lazy_t; ortho : Orthonormalization.t lazy_t; eN_ints : NucInt.t lazy_t; kin_ints : KinInt.t lazy_t; @@ -17,6 +18,7 @@ type t = let basis t = t.basis let overlap t = Lazy.force t.overlap +let multipole t = Lazy.force t.multipole let ortho t = Lazy.force t.ortho let eN_ints t = Lazy.force t.eN_ints let kin_ints t = Lazy.force t.kin_ints @@ -62,7 +64,11 @@ let make ~cartesian ~basis ?f12 nuclei = ScreenedERI.of_basis basis ) in - { basis ; overlap ; ortho ; eN_ints ; kin_ints ; ee_ints ; + let multipole = lazy ( + Multipole.of_basis basis + ) in + + { basis ; overlap ; multipole ; ortho ; eN_ints ; kin_ints ; ee_ints ; ee_lr_ints ; f12_ints ; f12_over_r12_ints ; cartesian ; } diff --git a/Basis/AOBasis.mli b/Basis/AOBasis.mli index 328bf4d..d00131b 100644 --- a/Basis/AOBasis.mli +++ b/Basis/AOBasis.mli @@ -10,6 +10,9 @@ val basis : t -> Basis.t val overlap : t -> Overlap.t (** Overlap matrix *) +val multipole : t -> Multipole.t +(** Multipole matrices *) + val ortho : t -> Orthonormalization.t (** Orthonormalization matrix of the overlap *) diff --git a/Basis/Basis.ml b/Basis/Basis.ml index e831d8f..9a3302f 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -93,6 +93,12 @@ let of_nuclei_and_basis_filename ~nuclei filename = in of_nuclei_and_general_basis nuclei general_basis +let of_nuclei_and_basis_string ~nuclei str = + let general_basis = + GeneralBasis.of_string str + in + of_nuclei_and_general_basis nuclei general_basis + let of_nuclei_and_basis_filenames ~nuclei filenames = let general_basis = diff --git a/Basis/Basis.mli b/Basis/Basis.mli index 83f010d..ea9297e 100644 --- a/Basis/Basis.mli +++ b/Basis/Basis.mli @@ -18,6 +18,11 @@ val of_nuclei_and_basis_filename : nuclei:Nuclei.t -> string -> t from a file. *) +val of_nuclei_and_basis_string : nuclei:Nuclei.t -> string -> t +(** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t} + from a string. + *) + val of_nuclei_and_basis_filenames : nuclei:Nuclei.t -> string list -> t (** Same as {!of_nuclei_and_general_basis}, but taking the {!GeneralBasis.t} from multiple files. diff --git a/Basis/GeneralBasis.ml b/Basis/GeneralBasis.ml index 4b8dcde..3f05319 100644 --- a/Basis/GeneralBasis.ml +++ b/Basis/GeneralBasis.ml @@ -17,23 +17,24 @@ exception No_shell exception Malformed_shell of string -let read_shell ic = +let read_shell line_stream = try let shell, n = try - let line = input_line ic in + let line = Stream.next line_stream in if String.trim line = "$END" then raise End_of_file; Scanf.sscanf line " %c %d " (fun shell n -> shell, n) with | End_of_file -> raise No_shell + | Stream.Failure -> raise No_shell | Scanf.Scan_failure m -> raise (Malformed_shell m) in let rec loop = function | 0 -> [] | i -> let contraction = - let line = (input_line ic) in + let line = Stream.next line_stream in try Scanf.sscanf line " %_d %f %f " (fun exponent coefficient -> { exponent ; coefficient }) with _ -> raise (Malformed_shell (Printf.sprintf @@ -50,36 +51,46 @@ let read_shell ic = -let read_element ic = +let read_element line_stream = try + let line = Stream.next line_stream in let element = - Scanf.sscanf (input_line ic) " %s " Element.of_string + Scanf.sscanf line " %s " Element.of_string in let rec loop () = - match read_shell ic with + match read_shell line_stream with | Some shell -> shell :: loop () | None -> [] in Some (element, Array.of_list (loop ()) ) with - | End_of_file -> None + | Stream.Failure -> None - -let read filename = - let ic = open_in filename in +let read_stream line_stream = let rec loop accu = try - match read_element ic with + match read_element line_stream with | Some e -> loop (e :: accu) | None -> accu with Element.ElementError _ -> loop accu in - loop [] + loop [] + +let read filename = + let ic = open_in filename in + let line_stream = + Stream.from (fun _ -> + try Some (input_line ic) + with End_of_file -> None ) + in + let result = read_stream line_stream in + close_in ic; + result let combine basis_list = @@ -124,3 +135,8 @@ let to_string (name, contracted_shell_array) = Printf.sprintf "%s\n%s" name (string_of_contracted_shell_array contracted_shell_array) +let of_string input_string = + String.split_on_char '\n' input_string + |> Stream.of_list + |> read_stream + diff --git a/Basis/GeneralBasis.mli b/Basis/GeneralBasis.mli index 138022a..2e0d540 100644 --- a/Basis/GeneralBasis.mli +++ b/Basis/GeneralBasis.mli @@ -46,7 +46,7 @@ exception Malformed_shell of string val read : string -> t -(** Reads a basis set file and return an association list where +(** Reads a basis-set file and return an association list where the key is an {!Element.t} and the value is the parsed basis set. *) @@ -65,8 +65,10 @@ val read_many : string list -> t -val read_element : in_channel -> element_basis option -(** Reads an element from the input channel [ic]. For example, +val read_element : string Stream.t -> element_basis option +(** Reads an element from the input [string Stream]. The [string Stream] is a +stream of lines, like a text file split on the end-of-line character. +For example, {[ HYDROGEN S 3 @@ -95,8 +97,10 @@ Some *) -val read_shell : in_channel -> general_contracted_shell option -(** Reads a shell from the input channel [ic]. For example, +val read_shell : string Stream.t -> general_contracted_shell option +(** Reads a shell from the input [string Stream]. The [string Stream] is a +stream of lines, like a text file split on the end-of-line character. +For example, {[ S 3 1 13.0100000 0.0196850 @@ -119,4 +123,6 @@ Some val to_string : string * (general_contracted_shell array) -> string (** Pretty-prints the basis set of an {Element.t}. *) +val of_string : string -> t +(** Reads a GAMESS-formatted string. *) diff --git a/Basis/Multipole.ml b/Basis/Multipole.ml new file mode 100644 index 0000000..b0289bf --- /dev/null +++ b/Basis/Multipole.ml @@ -0,0 +1,199 @@ +open Util +open Constants +open Lacaml.D + +type t = Mat.t array +(* +[| "x"; "y"; "z"; "x2"; "y2"; "z2" |] +*) + +module Am = AngularMomentum +module Bs = Basis +module Co = Coordinate +module Cs = ContractedShell +module Csp = ContractedShellPair +module Po = Powers +module Psp = PrimitiveShellPair + +let multiply a b = + let n = Mat.dim1 a in + let c = Mat.create n n in + Mat.cpab c a b; + c + +let x0 t = t.(0) +let y0 t = t.(1) +let z0 t = t.(2) +let x1 t = t.(3) +let y1 t = t.(4) +let z1 t = t.(5) +let x2 t = t.(6) +let y2 t = t.(7) +let z2 t = t.(8) + +let matrix_x t = multiply (x1 t) @@ multiply (y0 t) (z0 t) +let matrix_y t = multiply (x0 t) @@ multiply (y1 t) (z0 t) +let matrix_z t = multiply (x0 t) @@ multiply (y0 t) (z1 t) +let matrix_x2 t = multiply (x2 t) @@ multiply (y0 t) (z0 t) +let matrix_y2 t = multiply (x0 t) @@ multiply (y2 t) (z0 t) +let matrix_z2 t = multiply (x0 t) @@ multiply (y0 t) (z2 t) +let matrix_xy t = multiply (x1 t) @@ multiply (y1 t) (z0 t) +let matrix_yz t = multiply (x0 t) @@ multiply (y1 t) (z1 t) +let matrix_zx t = multiply (x1 t) @@ multiply (y0 t) (z1 t) + +let cutoff = integrals_cutoff + +let to_powers x = + let open Zkey in + match to_powers x with + | Six x -> x + | _ -> assert false + +(** Computes all the integrals of the contracted shell pair *) +let contracted_class shell_a shell_b : float Zmap.t array = + + match Csp.make shell_a shell_b with + | None -> Array.init 9 (fun _ -> Zmap.create 0) + | Some shell_p -> + begin + + (* Pre-computation of integral class indices *) + let class_indices = Csp.zkey_array shell_p in + + let contracted_class = + Array.init 9 (fun _ -> Array.make (Array.length class_indices) 0.) + in + + let a_minus_b = + Csp.a_minus_b shell_p + in + let norm_coef_scales = + Csp.norm_scales shell_p + in + + (* Compute all integrals in the shell for each pair of significant shell pairs *) + + let xyz_of_int k = + match k with + | 0 -> Co.X + | 1 -> Co.Y + | _ -> Co.Z + in + + List.iter (fun (coef_prod, psp) -> + (** Screening on the product of coefficients *) + if (abs_float coef_prod) > 1.e-6*.cutoff then + begin + let expo_inv = Psp.exponent_inv psp + and center_pa = Psp.center_minus_a psp + and xa = Co.get Co.X @@ Cs.center shell_a + and ya = Co.get Co.Y @@ Cs.center shell_a + and za = Co.get Co.Z @@ Cs.center shell_a + in + + Array.iteri (fun i key -> + let (angMomA, angMomB) = to_powers key in + (* 1D Overlap *) + let f k = + let xyz = xyz_of_int k in + Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB) + expo_inv + (Co.get xyz a_minus_b, Co.get xyz center_pa) + in + (* 1D *) + let g k = + let xyz = xyz_of_int k in + Overlap_primitives.hvrr (Po.get xyz angMomA + 1, Po.get xyz angMomB) + expo_inv + (Co.get xyz a_minus_b, Co.get xyz center_pa) + in + (* 1D *) + let h k = + let xyz = xyz_of_int k in + Overlap_primitives.hvrr (Po.get xyz angMomA + 2, Po.get xyz angMomB) + expo_inv + (Co.get xyz a_minus_b, Co.get xyz center_pa) + in + let norm = norm_coef_scales.(i) in + let f0, f1, f2, g0, g1, g2, h0, h1, h2 = + f 0, f 1, f 2, g 0, g 1, g 2, h 0, h 1, h 2 + in + let x = g0 +. f0 *. xa in + let y = g1 +. f1 *. ya in + let z = g2 +. f2 *. za in + let x2 = h0 +. xa *. (2. *. x -. xa *. f0) in + let y2 = h1 +. ya *. (2. *. y -. ya *. f1) in + let z2 = h2 +. za *. (2. *. z -. za *. f2) in + let c = contracted_class in + let d = coef_prod *. norm in + c.(0).(i) <- c.(0).(i) +. d *. f0 ; + c.(1).(i) <- c.(1).(i) +. d *. f1 ; + c.(2).(i) <- c.(2).(i) +. d *. f2 ; + c.(3).(i) <- c.(3).(i) +. d *. x ; + c.(4).(i) <- c.(4).(i) +. d *. y ; + c.(5).(i) <- c.(5).(i) +. d *. z ; + c.(6).(i) <- c.(6).(i) +. d *. x2 ; + c.(7).(i) <- c.(7).(i) +. d *. y2 ; + c.(8).(i) <- c.(8).(i) +. d *. z2 ; + ) class_indices + end + ) (Csp.coefs_and_shell_pairs shell_p); + let result = + Array.map (fun c -> Zmap.create (Array.length c) ) contracted_class + in + for j=0 to Array.length result do + let rj = result.(j) in + let cj = contracted_class.(j) in + Array.iteri (fun i key -> Zmap.add rj key cj.(i)) class_indices + done; + result + end + + +(** Create multipole matrices *) +let of_basis basis = + let to_powers x = + let open Zkey in + match to_powers x with + | Three x -> x + | _ -> assert false + in + + let n = Bs.size basis + and shell = Bs.contracted_shells basis + in + + let result = Array.init 9 (fun _ -> Mat.create n n) in + for j=0 to (Array.length shell) - 1 do + for i=0 to j do + (* Compute all the integrals of the class *) + let cls = + contracted_class shell.(i) shell.(j) + in + + for k=0 to 8 do + Array.iteri (fun j_c powers_j -> + let j_c = Cs.index shell.(j) + j_c + 1 in + let xj = to_powers powers_j in + Array.iteri (fun i_c powers_i -> + let i_c = Cs.index shell.(i) + i_c + 1 in + let xi = to_powers powers_i in + let key = + Zkey.of_powers_six xi xj + in + let value = + try Zmap.find cls.(k) key + with Not_found -> 0. + in + result.(k).{i_c,j_c} <- value; + result.(k).{j_c,i_c} <- value; + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(i)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom shell.(j)))) + done; + done; + done; + for k=0 to 8 do + Mat.detri result.(k); + done; + result + diff --git a/Basis/Multipole.mli b/Basis/Multipole.mli new file mode 100644 index 0000000..64d09de --- /dev/null +++ b/Basis/Multipole.mli @@ -0,0 +1,34 @@ +(** Multipole atomic integrals: + +{% $$ \langle \chi_i | x | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | y | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | z | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %} +{% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} + +*) + +open Lacaml.D + +type t + +val matrix_x : t -> Mat.t +(** {% $$ \langle \chi_i | x | \chi_j \rangle $$ %} *) + +val matrix_y : t -> Mat.t +(** {% $$ \langle \chi_i | y | \chi_j \rangle $$ %} *) + +val matrix_z : t -> Mat.t +(** {% $$ \langle \chi_i | z | \chi_j \rangle $$ %} *) + +val matrix_x2 : t -> Mat.t +(** {% $$ \langle \chi_i | x^2 | \chi_j \rangle $$ %} *) + +val matrix_y2 : t -> Mat.t +(** {% $$ \langle \chi_i | y^2 | \chi_j \rangle $$ %} *) + +val matrix_z2 : t -> Mat.t +(** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *) + +val of_basis : Basis.t -> t diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index f73dfa0..0eb858c 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -18,77 +18,77 @@ module Psp = PrimitiveShellPair let cutoff = integrals_cutoff -let to_powers x = +let to_powers x = let open Zkey in match to_powers x with | Six x -> x | _ -> assert false (** Computes all the overlap integrals of the contracted shell pair *) -let contracted_class shell_a shell_b : float Zmap.t = +let contracted_class shell_a shell_b : float Zmap.t = match Csp.make shell_a shell_b with - | Some shell_p -> + | None -> Zmap.create 0 + | Some shell_p -> begin - (* Pre-computation of integral class indices *) - let class_indices = Csp.zkey_array shell_p in + (* Pre-computation of integral class indices *) + let class_indices = Csp.zkey_array shell_p in - let contracted_class = - Array.make (Array.length class_indices) 0. - in + let contracted_class = + Array.make (Array.length class_indices) 0. + in - let a_minus_b = - Csp.a_minus_b shell_p - in - let norm_coef_scales = - Csp.norm_scales shell_p - in + let a_minus_b = + Csp.a_minus_b shell_p + in + let norm_coef_scales = + Csp.norm_scales shell_p + in - (* Compute all integrals in the shell for each pair of significant shell pairs *) + (* Compute all integrals in the shell for each pair of significant shell pairs *) - let xyz_of_int k = - match k with - | 0 -> Co.X - | 1 -> Co.Y - | _ -> Co.Z - in + let xyz_of_int k = + match k with + | 0 -> Co.X + | 1 -> Co.Y + | _ -> Co.Z + in - List.iter (fun (coef_prod, psp) -> - (** Screening on thr product of coefficients *) - if (abs_float coef_prod) > 1.e-3*.cutoff then - begin - let expo_inv = Psp.exponent_inv psp - and center_pa = Psp.center_minus_a psp - in - - Array.iteri (fun i key -> - let (angMomA,angMomB) = to_powers key in - let f k = - let xyz = xyz_of_int k in - Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB) - expo_inv - (Co.get xyz a_minus_b, - Co.get xyz center_pa) + List.iter (fun (coef_prod, psp) -> + (** Screening on the product of coefficients *) + if (abs_float coef_prod) > 1.e-6*.cutoff then + begin + let expo_inv = Psp.exponent_inv psp + and center_pa = Psp.center_minus_a psp in - let norm = norm_coef_scales.(i) 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 - ) (Csp.coefs_and_shell_pairs shell_p); - let result = - Zmap.create (Array.length contracted_class) - in - Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; - result - end - | None -> Zmap.create 0 + + Array.iteri (fun i key -> + let (angMomA,angMomB) = to_powers key in + let f k = + let xyz = xyz_of_int k in + Overlap_primitives.hvrr (Po.get xyz angMomA, Po.get xyz angMomB) + expo_inv + (Co.get xyz a_minus_b, + Co.get xyz center_pa) + in + let norm = norm_coef_scales.(i) 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 + ) (Csp.coefs_and_shell_pairs shell_p); + let result = + Zmap.create (Array.length contracted_class) + in + Array.iteri (fun i key -> Zmap.add result key contracted_class.(i)) class_indices; + result + end (** Create overlap matrix *) let of_basis basis = - let to_powers x = + let to_powers x = let open Zkey in match to_powers x with | Three x -> x @@ -104,7 +104,7 @@ let of_basis basis = for i=0 to j do (* Compute all the integrals of the class *) let cls = - contracted_class shell.(i) shell.(j) + contracted_class shell.(i) shell.(j) in Array.iteri (fun j_c powers_j -> @@ -113,11 +113,11 @@ let of_basis basis = Array.iteri (fun i_c powers_i -> let i_c = Cs.index shell.(i) + i_c + 1 in let xi = to_powers powers_i in - let key = + let key = Zkey.of_powers_six xi xj in - let value = - try Zmap.find cls key + let value = + try Zmap.find cls key with Not_found -> 0. in result.{i_c,j_c} <- value; @@ -132,7 +132,7 @@ let of_basis basis = (** Create mixed overlap matrix *) let of_basis_pair first_basis second_basis = - let to_powers x = + let to_powers x = let open Zkey in match to_powers x with | Three x -> x @@ -140,7 +140,7 @@ let of_basis_pair first_basis second_basis = in let n = Bs.size first_basis - and m = Bs.size second_basis + and m = Bs.size second_basis and first = Bs.contracted_shells first_basis and second = Bs.contracted_shells second_basis in @@ -150,7 +150,7 @@ let of_basis_pair first_basis second_basis = for i=0 to (Array.length first) - 1 do (* Compute all the integrals of the class *) let cls = - contracted_class first.(i) second.(j) + contracted_class first.(i) second.(j) in Array.iteri (fun j_c powers_j -> @@ -159,11 +159,11 @@ let of_basis_pair first_basis second_basis = Array.iteri (fun i_c powers_i -> let i_c = Cs.index first.(i) + i_c + 1 in let xi = to_powers powers_i in - let key = + let key = Zkey.of_powers_six xi xj in - let value = - try Zmap.find cls key + let value = + try Zmap.find cls key with Not_found -> 0. in result.{i_c,j_c} <- value; @@ -180,7 +180,7 @@ let to_file ~filename overlap = let oc = open_out filename in let n = - Mat.dim1 overlap + Mat.dim1 overlap in for j=1 to n do diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index 31a4ebd..33da7e4 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -3,11 +3,7 @@ open Xyz_ast type t = (Element.t * Coordinate.t) array -let of_xyz_file filename = - let lexbuf = - let ic = open_in filename in - Lexing.from_channel ic - in +let of_xyz_lexbuf lexbuf = let data = Xyz_parser.input Nuclei_lexer.read_all lexbuf in @@ -24,6 +20,30 @@ let of_xyz_file filename = |> Array.of_list +let of_xyz_string buffer = + Zmatrix.of_string buffer + |> Zmatrix.to_xyz + |> Array.map (fun (e,x,y,z) -> + (e, Coordinate.(angstrom_to_bohr @@ make_angstrom { x ; y ; z} )) + ) + + +let of_xyz_string input_string = + Lexing.from_string input_string + |> of_xyz_lexbuf + + +let of_xyz_file filename = + let ic = open_in filename in + let lexbuf = + Lexing.from_channel ic + in + let result = + of_xyz_lexbuf lexbuf + in + close_in ic; + result + let of_zmt_string buffer = Zmatrix.of_string buffer diff --git a/Nuclei/Nuclei.mli b/Nuclei/Nuclei.mli index 4d216cb..9944784 100644 --- a/Nuclei/Nuclei.mli +++ b/Nuclei/Nuclei.mli @@ -5,6 +5,9 @@ of tuples ({!Element.t}, {!Coordinate.t}). type t = (Element.t * Coordinate.t) array +val of_xyz_string : string -> t +(** Create from a string, in [xyz] format. *) + val of_xyz_file : string -> t (** Create from a file, in [xyz] format. *)