From ad9107ec8596fa8d3fbe4f97df39abaa83e4082f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 17 Sep 2014 23:47:13 +0200 Subject: [PATCH] Working on the basis set in xyz->ezfio --- ocaml/Basis.ml | 16 +++++++++------ ocaml/Basis.mli | 17 ++++++++++++++++ ocaml/Long_basis.ml | 10 ++++++---- ocaml/Long_basis.mli | 16 +++++++++++++++ ocaml/Makefile | 5 +++-- ocaml/qp_create_ezfio_from_xyz.ml | 33 +++++++++++++++++++++++++++++-- ocaml/qptypes_generator.ml | 7 +++++++ ocaml/test_basis.ml | 9 ++++++--- ocaml/test_gto.ml | 8 ++++---- 9 files changed, 100 insertions(+), 21 deletions(-) create mode 100644 ocaml/Basis.mli create mode 100644 ocaml/Long_basis.mli diff --git a/ocaml/Basis.ml b/ocaml/Basis.ml index 7fb4736d..ba575761 100644 --- a/ocaml/Basis.ml +++ b/ocaml/Basis.ml @@ -1,13 +1,14 @@ open Core.Std;; +open Qptypes;; -type t = Gto.t list;; +type t = (Gto.t * Atom_number.t) list;; (** Read all the basis functions of an element *) -let read in_channel = +let read in_channel at_number = let rec read result = try let gto = Gto.read_one in_channel in - read (gto::result) + read ( (gto,at_number)::result) with | Gto.End_Of_Basis -> List.rev result in read [] @@ -28,12 +29,15 @@ let find in_channel element = !element_read ;; -let read_element in_channel element = +(** Read an element from the file *) +let read_element in_channel at_number element = ignore (find in_channel element) ; - read in_channel ; + read in_channel at_number ; ;; let to_string b = - List.map ~f:Gto.to_string b + List.map ~f:(fun (g,n) -> + let n = Atom_number.to_int n in + (Int.to_string n)^":"^(Gto.to_string g)) b |> String.concat ~sep:"\n" ;; diff --git a/ocaml/Basis.mli b/ocaml/Basis.mli new file mode 100644 index 00000000..b9465d8c --- /dev/null +++ b/ocaml/Basis.mli @@ -0,0 +1,17 @@ +open Qptypes;; + +type t = (Gto.t * Atom_number.t) list + +(** Read all the basis functions of an element and set the number of the + * atom *) +val read : in_channel -> Atom_number.t -> (Gto.t * Atom_number.t) list + +(** Find an element in the basis set file *) +val find : in_channel -> Element.t -> Element.t + +(** Read the basis of an element from the file *) +val read_element : + in_channel -> Atom_number.t -> Element.t -> (Gto.t * Atom_number.t) list + +(** Convert the basis to a string *) +val to_string : (Gto.t * Atom_number.t) list -> string diff --git a/ocaml/Long_basis.ml b/ocaml/Long_basis.ml index b602a630..b5a600f9 100644 --- a/ocaml/Long_basis.ml +++ b/ocaml/Long_basis.ml @@ -1,15 +1,16 @@ open Core.Std;; +open Qptypes;; -type t = (Symmetry.Xyz.t * Gto.t) list;; +type t = (Symmetry.Xyz.t * Gto.t * Atom_number.t ) list;; let of_basis b = let rec do_work accu = function | [] -> accu - | g::tail -> + | (g,n)::tail -> begin let new_accu = Symmetry.Xyz.of_symmetry g.Gto.sym - |> List.map ~f:(fun x-> (x,g)) + |> List.map ~f:(fun x-> (x,g,n)) in do_work (new_accu@accu) tail end @@ -20,7 +21,8 @@ let of_basis b = let to_string b = - List.map ~f:(fun (x,y) -> + List.map ~f:(fun (x,y,z) -> + (Int.to_string (Atom_number.to_int z))^":"^ (Symmetry.Xyz.to_string x)^" "^(Gto.to_string y) ) b |> String.concat ~sep:"\n" diff --git a/ocaml/Long_basis.mli b/ocaml/Long_basis.mli new file mode 100644 index 00000000..934505eb --- /dev/null +++ b/ocaml/Long_basis.mli @@ -0,0 +1,16 @@ +open Qptypes;; + +(** A long basis is a basis set where + * all the P orbitals are converted to x, y, z + * all the D orbitals are converted to xx, xy, xz, yy, yx + * etc +*) +type t = (Symmetry.Xyz.t * Gto.t * Atom_number.t) list + +(** Transform a basis to a long basis *) +val of_basis : + (Gto.t * Atom_number.t) list -> (Symmetry.Xyz.t * Gto.t * Atom_number.t) list + +(** Convert the basis into its string representation *) +val to_string : + (Symmetry.Xyz.t * Gto.t * Atom_number.t) list -> string diff --git a/ocaml/Makefile b/ocaml/Makefile index 15f4baa9..ba3ff720 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -13,7 +13,7 @@ endif LIBS= PKGS= OCAMLCFLAGS=-g -OCAMLBUILD=ocamlbuild -cflags $(OCAMLCFLAGS) -lflags -g +OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags -g MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml MLIFILES=$(wildcard *.mli) ALL_TESTS=$(patsubst %.ml,%.byte,$(wildcard test_*.ml)) @@ -26,7 +26,8 @@ default: $(ALL_TESTS) $(ALL_EXE) executables: $(ALL_EXE) %.inferred.mli: $(MLFILES) - $(OCAMLBUILD) $*.inferred.mli -cflags -i -use-ocamlfind $(PKGS) + $(OCAMLBUILD) $*.inferred.mli -use-ocamlfind $(PKGS) + mv _build/$*.inferred.mli . %.byte: $(MLFILES) $(MLIFILES) $(OCAMLBUILD) $*.byte -use-ocamlfind $(PKGS) diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index a4716075..8918fbf8 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -81,8 +81,37 @@ File : %s\n" c m b xyz_file; ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords); (* Write Basis set *) - let basis = - let long_basis = Long_basis + let basis = + let rec do_work (accu:(Atom.t*Atom_number.t) list) (n:int) = function + | [] -> accu + | e::tail -> let new_accu = (e,(Atom_number.of_int n))::accu in + do_work new_accu (n+1) tail + in + do_work [] 1 nuclei + |> List.rev + |> List.map ~f:(fun (x,i) -> + Basis.read_element basis_channel i x.Atom.element) + |> List.concat + in + let long_basis = Long_basis.of_basis basis in + let ao_num = List.length long_basis in + Ezfio.set_ao_basis_ao_num ao_num; + let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc) + and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Atom_number.to_int n) + in + Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; + Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ; + +(* +ao_basis + ao_power integer (ao_basis_ao_num,3) + ao_prim_num_max integer = maxval(ao_basis_ao_prim_num) + ao_coef double precision (ao_basis_ao_num,ao_basis_ao_prim_num_max) + ao_expo double precision (ao_basis_ao_num,ao_basis_ao_prim_num_max) + +*) ;; diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index 0cab22d7..fb7ee875 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -28,6 +28,13 @@ let input_data = " * Non_empty_string : string assert (x <> \"\") ; +* Atom_number : int + assert (x > 0) ; + if (x > 1000) then + warning \"More than 1000 atoms\"; + if (Ezfio.has_nuclei_nucl_num ()) then + assert (x <= (Ezfio.get_nuclei_nucl_num ())); + * MO_number : int assert (x > 0) ; if (x > 1000) then diff --git a/ocaml/test_basis.ml b/ocaml/test_basis.ml index 278821f4..f8b056bf 100644 --- a/ocaml/test_basis.ml +++ b/ocaml/test_basis.ml @@ -1,5 +1,6 @@ open Core.Std;; open Qputils;; +open Qptypes;; let test_module () = @@ -13,12 +14,14 @@ let test_module () = Molecule.of_xyz_file xyz_file in + let nuclei = molecule.Molecule.nuclei in + let basis = - (Basis.read_element basis_channel Element.F) @ - (Basis.read_element basis_channel Element.F) + (Basis.read_element basis_channel (Atom_number.of_int 1) Element.F) @ + (Basis.read_element basis_channel (Atom_number.of_int 2) Element.F) in - Long_basis.of_basis basis + Long_basis.of_basis basis |> Long_basis.to_string |> print_endline ;; diff --git a/ocaml/test_gto.ml b/ocaml/test_gto.ml index bbfb8764..a39de97e 100644 --- a/ocaml/test_gto.ml +++ b/ocaml/test_gto.ml @@ -23,14 +23,14 @@ let test_gto_1 () = let test_gto_2 () = let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in ignore (input_line in_channel); - let basis = Basis.read in_channel in - List.iter basis ~f:(fun x-> Printf.printf "%s\n" (Gto.to_string x)) + let basis = Basis.read in_channel (Atom_number.of_int 1) in + List.iter basis ~f:(fun (x,n)-> Printf.printf "%d:%s\n" (Atom_number.to_int n) (Gto.to_string x)) ;; let test_gto () = let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in - let basis = Basis.read_element in_channel Element.C in - List.iter basis ~f:(fun x-> Printf.printf "%s\n" (Gto.to_string x)) + let basis = Basis.read_element in_channel (Atom_number.of_int 1) Element.C in + List.iter basis ~f:(fun (x,n)-> Printf.printf "%d:%s\n" (Atom_number.to_int n) (Gto.to_string x)) ;; let test_module () =