From 68b33e4b356a26723123e65ec90763627babb946 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 21 May 2021 16:42:48 +0200 Subject: [PATCH] Changed Symmetry into Angmom in OCaml --- ocaml/{Symmetry.ml => Angmom.ml} | 6 ++-- ocaml/{Symmetry.mli => Angmom.mli} | 0 ocaml/GaussianPrimitive.ml | 4 +-- ocaml/Gto.ml | 8 ++--- ocaml/Gto.mli | 2 +- ocaml/Input_ao_basis.ml | 18 ++++++------ ocaml/Long_basis.ml | 8 ++--- ocaml/Long_basis.mli | 8 ++--- ocaml/Primitive.mli | 4 +-- ocaml/qp_create_ezfio.ml | 7 +++-- src/basis/EZFIO.cfg | 47 ++++++++++++++++++++++++++++++ src/basis/NEED | 1 + src/basis/README.rst | 8 +++++ src/basis/basis.irp.f | 32 ++++++++++++++++++++ 14 files changed, 121 insertions(+), 32 deletions(-) rename ocaml/{Symmetry.ml => Angmom.ml} (95%) rename ocaml/{Symmetry.mli => Angmom.mli} (100%) create mode 100644 src/basis/EZFIO.cfg create mode 100644 src/basis/NEED create mode 100644 src/basis/README.rst create mode 100644 src/basis/basis.irp.f diff --git a/ocaml/Symmetry.ml b/ocaml/Angmom.ml similarity index 95% rename from ocaml/Symmetry.ml rename to ocaml/Angmom.ml index eb4b637b..ed13e8dc 100644 --- a/ocaml/Symmetry.ml +++ b/ocaml/Angmom.ml @@ -26,7 +26,7 @@ let of_string = function | "J" | "j" -> J | "K" | "k" -> K | "L" | "l" -> L - | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L, + | x -> raise (Failure ("Angmom should be S|P|D|F|G|H|I|J|K|L, not "^x^".")) let of_char = function @@ -40,7 +40,7 @@ let of_char = function | 'J' | 'j' -> J | 'K' | 'k' -> K | 'L' | 'l' -> L - | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L")) + | x -> raise (Failure ("Angmom should be S|P|D|F|G|H|I|J|K|L")) let to_l = function | S -> Positive_int.of_int 0 @@ -68,7 +68,7 @@ let of_l i = | 7 -> J | 8 -> K | 9 -> L - | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L")) + | x -> raise (Failure ("Angmom should be S|P|D|F|G|H|I|J|K|L")) type st = t diff --git a/ocaml/Symmetry.mli b/ocaml/Angmom.mli similarity index 100% rename from ocaml/Symmetry.mli rename to ocaml/Angmom.mli diff --git a/ocaml/GaussianPrimitive.ml b/ocaml/GaussianPrimitive.ml index d2144dec..c59a68de 100644 --- a/ocaml/GaussianPrimitive.ml +++ b/ocaml/GaussianPrimitive.ml @@ -2,14 +2,14 @@ open Qptypes open Sexplib.Std type t = - { sym : Symmetry.t ; + { sym : Angmom.t ; expo : AO_expo.t ; } [@@deriving sexp] let to_string p = let { sym = s ; expo = e } = p in Printf.sprintf "(%s, %22e)" - (Symmetry.to_string s) + (Angmom.to_string s) (AO_expo.to_float e) diff --git a/ocaml/Gto.ml b/ocaml/Gto.ml index 465e5627..200b4fa6 100644 --- a/ocaml/Gto.ml +++ b/ocaml/Gto.ml @@ -9,7 +9,7 @@ type fmt = | Gaussian type t = -{ sym : Symmetry.t ; +{ sym : Angmom.t ; lc : ((GaussianPrimitive.t * AO_coef.t) list) } [@@deriving sexp] @@ -47,7 +47,7 @@ let read_one in_channel = in let sym_str = String.sub buffer 0 2 in let n_str = String.sub buffer 2 ((String.length buffer)-2) in - let sym = Symmetry.of_string (String_ext.strip sym_str) in + let sym = Angmom.of_string (String_ext.strip sym_str) in let n = int_of_string (String_ext.strip n_str) in (* Read all the primitives *) let rec read_lines result = function @@ -82,7 +82,7 @@ let read_one in_channel = (** Write the GTO in Gamess format *) let to_string_gamess { sym = sym ; lc = lc } = let result = - Printf.sprintf "%s %3d" (Symmetry.to_string sym) (List.length lc) + Printf.sprintf "%s %3d" (Angmom.to_string sym) (List.length lc) in let rec do_work accu i = function | [] -> List.rev accu @@ -102,7 +102,7 @@ let to_string_gamess { sym = sym ; lc = lc } = (** Write the GTO in Gaussian format *) let to_string_gaussian { sym = sym ; lc = lc } = let result = - Printf.sprintf "%s %3d 1.00" (Symmetry.to_string sym) (List.length lc) + Printf.sprintf "%s %3d 1.00" (Angmom.to_string sym) (List.length lc) in let rec do_work accu i = function | [] -> List.rev accu diff --git a/ocaml/Gto.mli b/ocaml/Gto.mli index 91534ebe..ccbac7fe 100644 --- a/ocaml/Gto.mli +++ b/ocaml/Gto.mli @@ -5,7 +5,7 @@ type fmt = | Gaussian type t = - { sym : Symmetry.t ; + { sym : Angmom.t ; lc : (GaussianPrimitive.t * Qptypes.AO_coef.t) list; } [@@deriving sexp] diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index a73d641c..95d37a7a 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -9,7 +9,7 @@ module Ao_basis : sig ao_prim_num : AO_prim_number.t array; ao_prim_num_max : AO_prim_number.t; ao_nucl : Nucl_number.t array; - ao_power : Symmetry.Xyz.t array; + ao_power : Angmom.Xyz.t array; ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; @@ -32,7 +32,7 @@ end = struct ao_prim_num : AO_prim_number.t array; ao_prim_num_max : AO_prim_number.t; ao_nucl : Nucl_number.t array; - ao_power : Symmetry.Xyz.t array; + ao_power : Angmom.Xyz.t array; ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; @@ -87,7 +87,7 @@ end = struct if (data.(2*dim+i-1) > 0) then result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1)); done; - Array.map Symmetry.Xyz.of_string result + Array.map Angmom.Xyz.of_string result ;; let read_ao_coef () = @@ -133,7 +133,7 @@ end = struct let ao_num = AO_number.to_int b.ao_num in let gto_array = Array.init (AO_number.to_int b.ao_num) (fun i -> - let s = Symmetry.Xyz.to_symmetry b.ao_power.(i) in + let s = Angmom.Xyz.to_symmetry b.ao_power.(i) in let ao_prim_num = AO_prim_number.to_int b.ao_prim_num.(i) in let prims = List.init ao_prim_num (fun j -> let prim = { GaussianPrimitive.sym = s ; @@ -217,9 +217,9 @@ end = struct let ao_power = let l = Array.to_list ao_power in List.concat [ - (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.x) l) ; - (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.y) l) ; - (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.z) l) ] + (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.x) l) ; + (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.y) l) ; + (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.z) l) ] in Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; @@ -409,7 +409,7 @@ end = struct | [] -> [] | (i,n,x)::tail -> (Printf.sprintf " %5d %6d %-8s\n" i (Nucl_number.to_int n) - (Symmetry.Xyz.to_string x) + (Angmom.Xyz.to_string x) )::(do_work tail) in do_work l |> String.concat "" @@ -496,7 +496,7 @@ md5 = %s (b.ao_nucl |> Array.to_list |> list_map Nucl_number.to_string |> String.concat ", ") (b.ao_power |> Array.to_list |> list_map (fun x-> - "("^(Symmetry.Xyz.to_string x)^")" )|> String.concat ", ") + "("^(Angmom.Xyz.to_string x)^")" )|> String.concat ", ") (b.ao_coef |> Array.to_list |> list_map AO_coef.to_string |> String.concat ", ") (b.ao_expo |> Array.to_list |> list_map AO_expo.to_string diff --git a/ocaml/Long_basis.ml b/ocaml/Long_basis.ml index 6c88954a..a8ea3c66 100644 --- a/ocaml/Long_basis.ml +++ b/ocaml/Long_basis.ml @@ -2,7 +2,7 @@ open Qptypes open Qputils open Sexplib.Std -type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp] +type t = (Angmom.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp] let of_basis b = let rec do_work accu = function @@ -10,7 +10,7 @@ let of_basis b = | (g,n)::tail -> begin let new_accu = - Symmetry.Xyz.of_symmetry g.Gto.sym + Angmom.Xyz.of_symmetry g.Gto.sym |> List.rev_map (fun x-> (x,g,n)) in do_work (new_accu@accu) tail @@ -25,7 +25,7 @@ let to_basis b = | [] -> List.rev accu | (s,g,n)::tail -> let first_sym = - Symmetry.Xyz.of_symmetry g.Gto.sym + Angmom.Xyz.of_symmetry g.Gto.sym |> List.hd in let new_accu = @@ -42,7 +42,7 @@ let to_basis b = let to_string b = let middle = list_map (fun (x,y,z) -> "( "^((string_of_int (Nucl_number.to_int z)))^", "^ - (Symmetry.Xyz.to_string x)^", "^(Gto.to_string y) + (Angmom.Xyz.to_string x)^", "^(Gto.to_string y) ^" )" ) b |> String.concat ",\n" diff --git a/ocaml/Long_basis.mli b/ocaml/Long_basis.mli index 26009c10..2c41ad74 100644 --- a/ocaml/Long_basis.mli +++ b/ocaml/Long_basis.mli @@ -5,16 +5,16 @@ open Qptypes;; * all the D orbitals are converted to xx, xy, xz, yy, yx * etc *) -type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list [@@deriving sexp] +type t = (Angmom.Xyz.t * Gto.t * Nucl_number.t) list [@@deriving sexp] (** Transform a basis to a long basis *) val of_basis : - (Gto.t * Nucl_number.t) list -> (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list + (Gto.t * Nucl_number.t) list -> (Angmom.Xyz.t * Gto.t * Nucl_number.t) list (** Transform a long basis to a basis *) val to_basis : - (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> (Gto.t * Nucl_number.t) list + (Angmom.Xyz.t * Gto.t * Nucl_number.t) list -> (Gto.t * Nucl_number.t) list (** Convert the basis into its string representation *) val to_string : - (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> string + (Angmom.Xyz.t * Gto.t * Nucl_number.t) list -> string diff --git a/ocaml/Primitive.mli b/ocaml/Primitive.mli index f7d8809d..f63ecd8c 100644 --- a/ocaml/Primitive.mli +++ b/ocaml/Primitive.mli @@ -1,5 +1,5 @@ type t = -{ sym : Symmetry.t; +{ sym : Angmom.t; expo : Qptypes.AO_expo.t; } [@@deriving sexp] @@ -7,5 +7,5 @@ type t = val to_string : t -> string (** Creation *) -val of_sym_expo : Symmetry.t -> Qptypes.AO_expo.t -> t +val of_sym_expo : Angmom.t -> Qptypes.AO_expo.t -> t diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index f4c7f3b6..6a4a36a1 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -512,13 +512,14 @@ let run ?o b au c d m p cart xyz_file = let ao_num = List.length long_basis in Ezfio.set_ao_basis_ao_num ao_num; Ezfio.set_ao_basis_ao_basis b; + Ezfio.set_basis_basis b; let ao_prim_num = list_map (fun (_,g,_) -> List.length g.Gto.lc) long_basis and ao_nucl = list_map (fun (_,_,n) -> Nucl_number.to_int n) long_basis and ao_power= let l = list_map (fun (x,_,_) -> x) long_basis in - (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) l)@ - (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) l)@ - (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) l) + (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.x)) l)@ + (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.y)) l)@ + (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.z)) l) in let ao_prim_num_max = List.fold_left (fun s x -> if x > s then x diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg new file mode 100644 index 00000000..d560a291 --- /dev/null +++ b/src/basis/EZFIO.cfg @@ -0,0 +1,47 @@ +[basis] +type: character*(256) +doc: Name of the |AO| basis set +interface: ezfio + +[shell_num] +type: integer +doc: Number of shells +interface: ezfio, provider + +[shell_normalization_factor] +type: double precision +doc: Number of primitives per |AO| +size: (basis.shell_num) +interface: ezfio, provider + +[shell_prim_num] +type: integer +doc: Max number of primitives in a shell +size: (basis.shell_num) +interface: ezfio, provider + +[shell_prim_index] +type: integer +doc: Max number of primitives in a shell +size: (basis.shell_num) +interface: ezfio, provider + +[shell_nucl] +type: integer +doc: Index of the nucleus on which the shell is centered +size: (basis.shell_num) +interface: ezfio, provider + +[shell_coef] +type: double precision +doc: Primitive coefficients +size: (basis.shell_prim_num) +interface: ezfio, provider + +[shell_expo] +type: double precision +doc: Exponents in the shell +size: (basis.shell_prim_num) +interface: ezfio, provider + + diff --git a/src/basis/NEED b/src/basis/NEED new file mode 100644 index 00000000..d2066b18 --- /dev/null +++ b/src/basis/NEED @@ -0,0 +1 @@ +nuclei diff --git a/src/basis/README.rst b/src/basis/README.rst new file mode 100644 index 00000000..17493b4f --- /dev/null +++ b/src/basis/README.rst @@ -0,0 +1,8 @@ +====== +basis +====== + +This module contains the basis set information, which will then be used to build the atomic orbitals. + + + diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f new file mode 100644 index 00000000..2674b9f5 --- /dev/null +++ b/src/basis/basis.irp.f @@ -0,0 +1,32 @@ +BEGIN_PROVIDER [ double precision, basis_normalization_factor, (shell_num) ] + implicit none + BEGIN_DOC + ! Normalization factors of the shells + END_DOC + double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c + integer :: l, powA(3), nz + integer :: i,j,k + nz=100 + C_A(1) = 0.d0 + C_A(2) = 0.d0 + C_A(3) = 0.d0 + + do i=1,shell_num + powA(1) = shell_ang_mom(i) + powA(2) = 0 + powA(3) = 0 + + ! Normalization of the contracted basis functions + norm = 0.d0 + do j=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 + do k=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 + call overlap_gaussian_xyz(C_A,C_A,shell_prim_expo(j),shell_prim_expo(k), & + powA,powA,overlap_x,overlap_y,overlap_z,c,nz) + norm = norm+c*shell_prim_coef(j)*shell_prim_coef(k) + enddo + enddo + basis_normalization_factor(i) = basis_normalization_factor(i) * sqrt(norm) + + enddo + +END_PROVIDER