From 0fb8c3b7ae614d7f16463b71536144a5d5950195 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 18 Sep 2014 17:01:43 +0200 Subject: [PATCH] xyz to ezfio works --- ocaml/Atom.ml | 7 ++-- ocaml/Gto.ml | 6 ++-- ocaml/Makefile | 2 +- ocaml/Molecule.ml | 3 +- ocaml/Point3d.ml | 20 ++++++++--- ocaml/Qputils.ml | 10 ++++++ ocaml/Symmetry.ml | 5 +-- ocaml/Symmetry.mli | 18 ++++++++++ ocaml/qp_create_ezfio_from_xyz.ml | 58 +++++++++++++++++++++++++------ ocaml/qptypes_generator.ml | 18 ++++++---- ocaml/test_atom.ml | 2 +- ocaml/test_point3d.ml | 12 ++++--- robots.txt | 2 -- src/MOs/mos.irp.f | 35 ++++++++++++------- src/MOs/utils.irp.f | 1 + 15 files changed, 149 insertions(+), 50 deletions(-) create mode 100644 ocaml/Symmetry.mli delete mode 100644 robots.txt diff --git a/ocaml/Atom.ml b/ocaml/Atom.ml index b9ab2866..d7ff0fdb 100644 --- a/ocaml/Atom.ml +++ b/ocaml/Atom.ml @@ -26,7 +26,8 @@ type t = coord : Point3d.t ; } -let of_string s = +(** Read xyz coordinates of the atom with unit u *) +let of_string u s = let buffer = s |> String.split ~on:' ' |> List.filter ~f:(fun x -> x <> "") @@ -35,13 +36,13 @@ let of_string s = | [ name; charge; x; y; z ] -> { element = Element.of_string name ; charge = Charge.of_string charge ; - coord = Point3d.of_string (String.concat [x; y; z] ~sep:" ") + coord = Point3d.of_string u (String.concat [x; y; z] ~sep:" ") } | [ name; x; y; z ] -> let e = Element.of_string name in { element = e ; charge = Charge.of_int (Element.charge e); - coord = Point3d.of_string (String.concat [x; y; z] ~sep:" ") + coord = Point3d.of_string u (String.concat [x; y; z] ~sep:" ") } | _ -> raise (AtomError s) ;; diff --git a/ocaml/Gto.ml b/ocaml/Gto.ml index 88f5d874..8abc3837 100644 --- a/ocaml/Gto.ml +++ b/ocaml/Gto.ml @@ -6,7 +6,7 @@ exception End_Of_Basis type t = { sym : Symmetry.t ; - lc : ((Primitive.t * float) list) + lc : ((Primitive.t * AO_coef.t) list) } ;; @@ -56,7 +56,7 @@ let read_one in_channel = Primitive.expo = Positive_float.of_float (Float.of_string expo) } - and c = Float.of_string coef in + and c = AO_coef.of_float (Float.of_string coef) in read_lines ( (p,c)::result) (i-1) end | _ -> raise (GTO_Read_Failure line_buffer) @@ -68,7 +68,7 @@ let read_one in_channel = (** Transform the gto to a string *) let to_string { sym = sym ; lc = lc } = - let f (p,c) = Printf.sprintf "( %s, %f )" (Primitive.to_string p) c + let f (p,c) = Printf.sprintf "( %s, %f )" (Primitive.to_string p) (AO_coef.to_float c) in Printf.sprintf "[ %s : %s ]" (Symmetry.to_string sym) (String.concat (List.map ~f:f lc) ~sep:", ") diff --git a/ocaml/Makefile b/ocaml/Makefile index ba3ff720..80f19054 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -12,7 +12,7 @@ endif LIBS= PKGS= -OCAMLCFLAGS=-g +OCAMLCFLAGS=-g OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags -g MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml MLIFILES=$(wildcard *.mli) diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index df393168..7e1cde47 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -76,10 +76,11 @@ let to_string m = let of_xyz_string ?(charge=0) ?(multiplicity=(Multiplicity.of_int 1)) + ?(units=Point3d.Angstrom) s = let l = String.split s ~on:'\n' |> List.filter ~f:(fun x -> x <> "") - |> List.map ~f:Atom.of_string + |> List.map ~f:(fun x -> Atom.of_string units x) in let ne = ( get_charge { nuclei=l ; diff --git a/ocaml/Point3d.ml b/ocaml/Point3d.ml index 58b39d45..ecc4963b 100644 --- a/ocaml/Point3d.ml +++ b/ocaml/Point3d.ml @@ -1,21 +1,33 @@ open Core.Std;; +type units = +| Bohr +| Angstrom +;; + type t = { x : float ; y : float ; z : float ; } -let of_string s = +(** Read x y z coordinates in string s with units u *) +let of_string u s = + let f = + begin match u with + | Bohr -> 1. + | Angstrom -> 1. /. 0.52917721092 + end + in let l = s |> String.split ~on:' ' |> List.filter ~f:(fun x -> x <> "") |> List.map ~f:Float.of_string |> Array.of_list in - { x = l.(0) ; - y = l.(1) ; - z = l.(2) } + { x = l.(0) *. f ; + y = l.(1) *. f ; + z = l.(2) *. f } ;; diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index ca042fc2..12cb7e06 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -1,4 +1,14 @@ let (/) (a:string) (b:string) = a^"/"^b;; +let rec transpose = function +| [] -> [] +| []::tail -> transpose tail +| (x::t1)::t2 -> + let new_head = (x::(List.map List.hd t2)) + and new_tail = (transpose (t1 :: (List.map List.tl t2) )) + in + new_head @ new_tail +;; + diff --git a/ocaml/Symmetry.ml b/ocaml/Symmetry.ml index aac3ed18..828d0dfa 100644 --- a/ocaml/Symmetry.ml +++ b/ocaml/Symmetry.ml @@ -58,7 +58,9 @@ type st = t ;; module Xyz : sig - type t + type t = { x: Positive_int.t ; + y: Positive_int.t ; + z: Positive_int.t } val of_string : string -> t val to_string : t -> string val get_l : t -> Positive_int.t @@ -67,7 +69,6 @@ end = struct type t = { x: Positive_int.t ; y: Positive_int.t ; z: Positive_int.t } - type state_type = Null | X | Y | Z (** Builds an XYZ triplet from a string. diff --git a/ocaml/Symmetry.mli b/ocaml/Symmetry.mli new file mode 100644 index 00000000..086a8a1a --- /dev/null +++ b/ocaml/Symmetry.mli @@ -0,0 +1,18 @@ +open Qptypes;; + +type t = S | P | D | F | G | H | I | J | K | L +val to_string : t -> string +val of_string : string -> t +val of_char : char -> t +val to_l : t -> Positive_int.t +type st = t +module Xyz : + sig + type t = { x: Positive_int.t ; + y: Positive_int.t ; + z: Positive_int.t } + val of_string : string -> t + val to_string : t -> string + val get_l : t -> Positive_int.t + val of_symmetry : st -> t list + end diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index 8918fbf8..cee24011 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -98,21 +98,59 @@ File : %s\n" c m b xyz_file; 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) + and ao_power= + let l = List.map long_basis ~f:(fun (x,_,_) -> x) in + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@ + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@ + (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) ) + in + let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x -> + if x > s then x + else s) ao_prim_num + in + let gtos = List.map long_basis ~f:(fun (_,x,_) -> x) in + + let create_expo_coef ec = + let coefs = + begin match ec with + | `Coefs -> List.map gtos ~f:(fun x-> + List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) ) + | `Expos -> List.map gtos ~f:(fun x-> + List.map x.Gto.lc ~f:(fun (prim,_) -> Positive_float.to_float + prim.Primitive.expo) ) + end + in + let rec get_n n accu = function + | [] -> List.rev accu + | h::tail -> + let y = + begin match List.nth h n with + | Some x -> x + | None -> 0. + end + in + get_n n (y::accu) tail + in + let rec build accu = function + | n when n=ao_prim_num_max -> accu + | n -> build ( accu @ (get_n n [] coefs) ) (n+1) + in + build [] 0 + in + + let ao_coef = create_expo_coef `Coefs + and ao_expo = create_expo_coef `Expos 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) - -*) - + Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; + Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ; + Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list + ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ; ;; let command = diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index fb7ee875..2c3f557a 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -54,7 +54,7 @@ let input_data = " if (x > 100) then warning \"N_int > 100\"; if (Ezfio.has_determinants_n_int ()) then - assert (x == (Ezfio.get_determinants_n_int ())); + assert (x = (Ezfio.get_determinants_n_int ())); * Det_number : int assert (x > 0) ; @@ -75,6 +75,10 @@ let input_data = " | _ -> raise (Failure \"Bit_kind should be (1|2|4|8).\") end; +* MO_coef : float + +* AO_coef : float + " ;; @@ -84,17 +88,19 @@ module %s : sig type t val to_%s : t -> %s val of_%s : %s -> t + val to_string : %s -> string end = struct type t = %s let to_%s x = x let of_%s x = ( %s x ) + let to_string x = %s.to_string x end " ;; let parse_input input= - print_string "let warning = print_string;;\n" ; + print_string "open Core.Std;;\nlet warning = print_string;;\n" ; let rec parse result = function | [] -> result | ( "" , "" )::tail -> parse result tail @@ -102,10 +108,10 @@ let parse_input input= let name , typ = String.lsplit2_exn ~on:':' t in let typ = String.strip typ - and name = String.strip name - in - let newstring = Printf.sprintf template name typ typ typ typ typ typ typ - ( String.strip text ) + and name = String.strip name in + let typ_cap = String.capitalize typ in + let newstring = Printf.sprintf template name typ typ typ typ typ typ typ typ + ( String.strip text ) typ_cap in List.rev (parse (newstring::result) tail ) in diff --git a/ocaml/test_atom.ml b/ocaml/test_atom.ml index 8e6c019b..02106e3f 100644 --- a/ocaml/test_atom.ml +++ b/ocaml/test_atom.ml @@ -1,5 +1,5 @@ let test_atom = let line = "C 6.0 1. 2. 3." in - let atom = Atom.of_string line in + let atom = Atom.of_string Point3d.Bohr line in print_string (Atom.to_string atom) ;; diff --git a/ocaml/test_point3d.ml b/ocaml/test_point3d.ml index 86c25c69..40e97ea9 100644 --- a/ocaml/test_point3d.ml +++ b/ocaml/test_point3d.ml @@ -1,13 +1,15 @@ let test_point3d_1 () = let input = "7.4950000 -0.1499810 0.5085570" in - let p3d = Point3d.of_string input in - print_string (Point3d.to_string p3d) + let p3d = Point3d.of_string Point3d.Angstrom input in + print_endline(Point3d.to_string p3d); + let p3d = Point3d.of_string Point3d.Bohr input in + print_endline (Point3d.to_string p3d) ;; let test_point3d () = - let p1 = Point3d.of_string "1. 2. 3." - and p2 = Point3d.of_string "-2. 1. 1.5" in + let p1 = Point3d.of_string Point3d.Bohr "1. 2. 3." + and p2 = Point3d.of_string Point3d.Bohr "-2. 1. 1.5" in Printf.printf "%f\n" (Point3d.distance p1 p2) ;; -test_point3d (); +test_point3d_1 (); diff --git a/robots.txt b/robots.txt deleted file mode 100644 index c2a49f4f..00000000 --- a/robots.txt +++ /dev/null @@ -1,2 +0,0 @@ -User-agent: * -Allow: / diff --git a/src/MOs/mos.irp.f b/src/MOs/mos.irp.f index 738d2c2a..1a0cea77 100644 --- a/src/MOs/mos.irp.f +++ b/src/MOs/mos.irp.f @@ -4,7 +4,13 @@ BEGIN_PROVIDER [ integer, mo_tot_num ] ! Total number of molecular orbitals and the size of the keys corresponding END_DOC PROVIDE ezfio_filename - call ezfio_get_mo_basis_mo_tot_num(mo_tot_num) + logical :: exists + call ezfio_has_mo_basis_mo_tot_num(exists) + if (exists) then + call ezfio_get_mo_basis_mo_tot_num(mo_tot_num) + else + mo_tot_num = ao_num + endif ASSERT (mo_tot_num > 0) END_PROVIDER @@ -42,18 +48,23 @@ END_PROVIDER endif ! Coefs - allocate(buffer(ao_num,mo_tot_num)) - buffer = 0.d0 - call ezfio_get_mo_basis_mo_coef(buffer) - do i=1,mo_tot_num - do j=1,ao_num - mo_coef(j,i) = buffer(j,i) + call ezfio_has_mo_basis_mo_coef(exists) + if (exists) then + allocate(buffer(ao_num,mo_tot_num)) + buffer = 0.d0 + call ezfio_get_mo_basis_mo_coef(buffer) + do i=1,mo_tot_num + do j=1,ao_num + mo_coef(j,i) = buffer(j,i) + enddo + do j=ao_num+1,ao_num_align + mo_coef(j,i) = 0.d0 + enddo enddo - do j=ao_num+1,ao_num_align - mo_coef(j,i) = 0.d0 - enddo - enddo - deallocate(buffer) + deallocate(buffer) + else + mo_coef = 0.d0 + endif END_PROVIDER diff --git a/src/MOs/utils.irp.f b/src/MOs/utils.irp.f index 5e7984ff..7e23f744 100644 --- a/src/MOs/utils.irp.f +++ b/src/MOs/utils.irp.f @@ -5,6 +5,7 @@ subroutine save_mos call system('$QPACKAGE_ROOT/scripts/save_current_mos.sh '//trim(ezfio_filename)) + call ezfio_set_mo_basis_mo_tot_num(mo_tot_num) call ezfio_set_mo_basis_mo_label(mo_label) allocate ( buffer(ao_num,mo_tot_num) ) buffer = 0.d0