Added qp_print.ml

This commit is contained in:
Anthony Scemama 2014-10-07 19:33:11 +02:00
parent 5d72e3d4ff
commit e7b81a11de
14 changed files with 282 additions and 60 deletions

View File

@ -2,24 +2,6 @@ open Core.Std;;
exception AtomError of string
module Charge : sig
type t
val to_float : t -> float
val to_int : t -> int
val to_string: t -> string
val of_float : float -> t
val of_int : int -> t
val of_string: string -> t
end = struct
type t = float
let to_float x = x
let to_int x = Float.to_int x
let to_string x = Float.to_string (to_float x)
let of_float x = x
let of_int i = Float.of_int i
let of_string s = Float.of_string s
end
type t =
{ element : Element.t ;
charge : Charge.t ;
@ -41,16 +23,16 @@ let of_string u s =
| [ name; x; y; z ] ->
let e = Element.of_string name in
{ element = e ;
charge = Charge.of_int (Element.charge e);
charge = Charge.of_int (Element.to_charge e);
coord = Point3d.of_string u (String.concat [x; y; z] ~sep:" ")
}
| _ -> raise (AtomError s)
;;
let to_string a =
let to_string u a =
[ Element.to_string a.element ;
Charge.to_string a.charge ;
Point3d.to_string a.coord ]
Point3d.to_string u a.coord ]
|> String.concat ?sep:(Some " ")
;;

17
ocaml/Charge.ml Normal file
View File

@ -0,0 +1,17 @@
open Core.Std;;
type t = float
let of_float x = x
let of_int i = Float.of_int i
let of_string s = Float.of_string s
let to_float x = x
let to_int x = Float.to_int x
let to_string x =
if x >= 0. then
Printf.sprintf "+%f" x
else
Printf.sprintf "%f" x
;;

8
ocaml/Charge.mli Normal file
View File

@ -0,0 +1,8 @@
type t = float
val to_float : t -> float
val to_int : t -> int
val to_string: t -> string
val of_float : float -> t
val of_int : int -> t
val of_string: string -> t

View File

@ -46,7 +46,7 @@ let of_string = function
| "Se" | "Selenium" -> Se
| "Br" | "Bromine" -> Br
| "Kr" | "Krypton" -> Kr
| x -> raise (ElementError ("Atom "^x^" unknown"))
| x -> raise (ElementError ("Element "^x^" unknown"))
;;
let to_string = function
@ -129,7 +129,7 @@ let to_long_string = function
| Kr -> "Krypton"
;;
let charge = function
let to_charge = function
| X -> 0
| H -> 1
| He -> 2
@ -169,3 +169,44 @@ let charge = function
| Kr -> 36
;;
let of_charge = function
| 0 -> X
| 1 -> H
| 2 -> He
| 3 -> Li
| 4 -> Be
| 5 -> B
| 6 -> C
| 7 -> N
| 8 -> O
| 9 -> F
| 10 -> Ne
| 11 -> Na
| 12 -> Mg
| 13 -> Al
| 14 -> Si
| 15 -> P
| 16 -> S
| 17 -> Cl
| 18 -> Ar
| 19 -> K
| 20 -> Ca
| 21 -> Sc
| 22 -> Ti
| 23 -> V
| 24 -> Cr
| 25 -> Mn
| 26 -> Fe
| 27 -> Co
| 28 -> Ni
| 29 -> Cu
| 30 -> Zn
| 31 -> Ga
| 32 -> Ge
| 33 -> As
| 34 -> Se
| 35 -> Br
| 36 -> Kr
| x -> raise (ElementError ("Element of charge "^(string_of_int x)^" unknown"))
;;

View File

@ -32,7 +32,7 @@ executables: $(ALL_EXE)
%.byte: $(MLFILES) $(MLIFILES)
rm -f -- $*
$(OCAMLBUILD) $*.byte -use-ocamlfind $(PKGS)
ln -s $*.native $*
ln -s $*.byte $*
%.native: $(MLFILES) $(MLIFILES)
rm -f -- $*

View File

@ -12,10 +12,10 @@ type t = {
let get_charge { nuclei ; elec_alpha ; elec_beta } =
let result = Positive_int.(to_int elec_alpha + to_int elec_beta) in
let rec nucl_charge = function
| a::rest -> Atom.(Charge.to_float a.charge) +. nucl_charge rest
| a::rest -> (Charge.to_float a.Atom.charge) +. nucl_charge rest
| [] -> 0.
in
nucl_charge nuclei -. (Float.of_int result)
Charge.of_float (nucl_charge nuclei -. (Float.of_int result))
;;
let get_multiplicity m =
@ -27,7 +27,7 @@ let get_nucl_num m =
;;
let name m =
let cm = Float.to_int (get_charge m) in
let cm = Charge.to_int (get_charge m) in
let c =
match cm with
| 0 -> ""
@ -70,13 +70,14 @@ let to_string m =
let { nuclei ; elec_alpha ; elec_beta } = m in
let n = List.length nuclei in
let title = name m in
[ Int.to_string n ; title ] @ (List.map ~f:Atom.to_string nuclei)
[ Int.to_string n ; title ] @ (List.map ~f:(fun x -> Atom.to_string
Units.Angstrom x) nuclei)
|> String.concat ~sep:"\n"
;;
let of_xyz_string
?(charge=0) ?(multiplicity=(Multiplicity.of_int 1))
?(units=Point3d.Angstrom)
?(units=Units.Angstrom)
s =
let l = String.split s ~on:'\n'
|> List.filter ~f:(fun x -> x <> "")
@ -86,8 +87,8 @@ let of_xyz_string
nuclei=l ;
elec_alpha=(Positive_int.of_int 0) ;
elec_beta=(Positive_int.of_int 0) }
|> Float.to_int
)- charge
|> Charge.to_int
) - charge
|> Positive_int.of_int
in
let (na,nb) = Multiplicity.to_alpha_beta ne multiplicity in

View File

@ -1,10 +1,5 @@
open Core.Std;;
type units =
| Bohr
| Angstrom
;;
type t = {
x : float ;
y : float ;
@ -13,11 +8,9 @@ type t = {
(** 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
let f = match u with
| Units.Bohr -> 1.
| Units.Angstrom -> Units.angstrom_to_bohr
in
let l = s
|> String.split ~on:' '
@ -40,8 +33,12 @@ let distance2 p1 p2 =
let distance p1 p2 = sqrt (distance2 p1 p2)
;;
let to_string p =
let to_string u p =
let f = match u with
| Units.Bohr -> 1.
| Units.Angstrom -> Units.bohr_to_angstrom
in
let { x=x ; y=y ; z=z } = p in
Printf.sprintf "%f %f %f" x y z
Printf.sprintf "%f %f %f" (x*.f) (y*.f) (z*.f)
;;

11
ocaml/Units.ml Normal file
View File

@ -0,0 +1,11 @@
type units =
| Bohr
| Angstrom
;;
let angstrom_to_bohr = 1. /. 0.52917721092
let bohr_to_angstrom = 0.52917721092
;;

View File

@ -18,14 +18,6 @@ let spec =
let run ?o b c m xyz_file =
(*
(* DEBUG *)
Printf.printf "Charge : %d
Multiplicity : %d
Basis : %s
File : %s\n" c m b xyz_file;
*)
(* Open basis set channel *)
let basis_channel =
In_channel.create

168
ocaml/qp_print.ml Normal file
View File

@ -0,0 +1,168 @@
open Qputils;;
open Qptypes;;
open Core.Std;;
type input_action =
| Basis
| Nuclei
| Charge
| Multiplicity
| Electrons
;;
let create_i_action = function
| "basis" -> Basis
| "nucl" -> Nuclei
| "charge" -> Charge
| "mult" -> Multiplicity
| "elec" -> Electrons
| _ -> raise (Failure "Action should be:
* basis
* nucl
* charge
* mult
* elec
")
;;
let spec =
let open Command.Spec in
empty
+> flag "i" (optional string)
~doc:"Prints input data"
+> flag "o" (optional string)
~doc:"Prints output data"
+> anon ("ezfio_file" %: string)
;;
let run_i ~action ezfio_filename =
let action = create_i_action action in
if (not (Sys.file_exists_exn ezfio_filename)) then
failwith (ezfio_filename^" does not exists");
Ezfio.set_file ezfio_filename;
let print_basis () =
Printf.printf "%s\n" (Ezfio.get_ao_basis_ao_basis ())
in
let compute_charge () =
let nucl_charge = Ezfio.((get_nuclei_nucl_charge ()).data)
|> Ezfio.flattened_ezfio_data |> Array.map ~f:(Float.to_int)
and n_alpha = Ezfio.get_electrons_elec_alpha_num ()
and n_beta = Ezfio.get_electrons_elec_beta_num ()
in Array.fold ~init:(-n_alpha-n_beta) ~f:(fun x y -> x+y) nucl_charge
in
let compute_multiplicity () =
let n_alpha = Positive_int.of_int (Ezfio.get_electrons_elec_alpha_num ())
and n_beta = Positive_int.of_int (Ezfio.get_electrons_elec_beta_num ())
in Multiplicity.of_alpha_beta n_alpha n_beta
in
let create_molecule () =
let nucl_num = Ezfio.get_nuclei_nucl_num ()
and nucl_charge = Ezfio.((get_nuclei_nucl_charge ()).data)
|> Ezfio.flattened_ezfio_data
and nucl_coord = Ezfio.((get_nuclei_nucl_coord ()).data )
|> Ezfio.flattened_ezfio_data
in
let nucl_label =
if (Ezfio.has_nuclei_nucl_label ()) then
Ezfio.((get_nuclei_nucl_label ()).data) |> Ezfio.flattened_ezfio_data
else
Array.map ~f:(fun x-> x
|> Float.to_int
|> Element.of_charge
|> Element.to_string ) nucl_charge
in
let buffer = ref "" in
for i=0 to (nucl_num-1) do
buffer := !buffer ^ (Printf.sprintf "%s %f %f %f %f\n"
nucl_label.(i)
nucl_charge.(i)
nucl_coord.(i)
nucl_coord.(i+nucl_num)
nucl_coord.(i+nucl_num+nucl_num)
)
done ;
let charge = compute_charge () in
let mult = compute_multiplicity () in
Molecule.of_xyz_string ~charge:charge ~multiplicity:mult !buffer
in
let print_nuclei () =
let molecule = create_molecule () in
print_endline (Molecule.to_string molecule)
and print_charge () =
let molecule = create_molecule () in
Printf.printf "%s" (Charge.to_string (Molecule.get_charge molecule))
and print_multiplicity () =
let molecule = create_molecule () in
Printf.printf "%s" (Multiplicity.to_string (Molecule.get_multiplicity
molecule))
and print_electrons () = ()
in
match action with
| Basis -> print_basis ()
| Nuclei -> print_nuclei ()
| Charge -> print_charge ()
| Multiplicity -> print_multiplicity ()
| Electrons -> print_electrons ()
;;
let run_o ~action ezfio_filename =
if (not (Sys.file_exists_exn ezfio_filename)) then
failwith (ezfio_filename^" does not exists");
(* Open EZFIO *)
Ezfio.set_file ezfio_filename;
;;
let command =
Command.basic
~summary: "Quantum Package command"
~readme:(fun () ->
"
Prints data contained into the EZFIO file.
Input data :
* basis
* nuclei
* charge
* multiplicity
* electrons
Output data :
*
")
spec
(fun i o ezfio_file () ->
match (i,o) with
| (Some i, None) -> run_i ~action:i ezfio_file
| (None, Some o) -> run_o ~action:o ezfio_file
| (Some _, Some _) ->
raise (Failure "Error : please specify -i or -o but not both.")
| (None, None) ->
raise (Failure "Error : please specify -i or -o.")
)
;;
let () =
Command.run command
;;

View File

@ -1,5 +1,5 @@
let test_atom =
let line = "C 6.0 1. 2. 3." in
let atom = Atom.of_string Point3d.Bohr line in
print_string (Atom.to_string atom)
let atom = Atom.of_string Units.Bohr line in
print_string (Atom.to_string Units.Angstrom atom)
;;

View File

@ -1,6 +1,6 @@
let test_module () =
let atom = Element.of_string "Cobalt" in
Printf.printf "%s %d\n" (Element.to_string atom) (Element.charge atom)
Printf.printf "%s %d\n" (Element.to_string atom) (Element.to_charge atom)
;;
test_module ();;

View File

@ -1,15 +1,19 @@
let test_point3d_1 () =
let input = "7.4950000 -0.1499810 0.5085570" in
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 p3d = Point3d.of_string Units.Angstrom input in
print_endline(Point3d.to_string Units.Angstrom p3d);
let p3d = Point3d.of_string Units.Bohr input in
print_endline (Point3d.to_string Units.Bohr p3d)
;;
let test_point3d () =
let p1 = Point3d.of_string Point3d.Bohr "1. 2. 3."
and p2 = Point3d.of_string Point3d.Bohr "-2. 1. 1.5" in
let p1 = Point3d.of_string Units.Bohr "1. 2. 3."
and p2 = Point3d.of_string Units.Bohr "-2. 1. 1.5" in
Printf.printf "%f\n" (Point3d.distance p1 p2);
let p1 = Point3d.of_string Units.Angstrom "1. 2. 3."
and p2 = Point3d.of_string Units.Angstrom "-2. 1. 1.5" in
Printf.printf "%f\n" (Point3d.distance p1 p2)
;;
test_point3d_1 ();
test_point3d ();

View File

@ -1,4 +1,5 @@
ao_basis
ao_basis character*(256)
ao_num integer
ao_prim_num integer (ao_basis_ao_num)
ao_nucl integer (ao_basis_ao_num)