From 9bffc0883bdc44a833a170f7a0491eedaa6747f9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 24 Aug 2014 20:00:26 +0200 Subject: [PATCH] Added molecule in ocaml --- ocaml/.gitignore | 5 ++ ocaml/Makefile | 26 +++++++++-- ocaml/README.rst | 6 +++ ocaml/atom.ml | 8 +++- ocaml/molecule.ml | 101 +++++++++++++++++++++++++++++++++++++++++ ocaml/multiplicity.ml | 35 ++++++++++++++ ocaml/point3d.ml | 2 +- ocaml/qptypes.ml | 8 ++-- ocaml/test_elements.ml | 4 +- ocaml/test_molecule.ml | 35 ++++++++++++++ src/Makefile.common | 1 + src/NEEDED_MODULES | 2 - 12 files changed, 219 insertions(+), 14 deletions(-) create mode 100644 ocaml/.gitignore create mode 100644 ocaml/README.rst create mode 100644 ocaml/molecule.ml create mode 100644 ocaml/multiplicity.ml create mode 100644 ocaml/test_molecule.ml diff --git a/ocaml/.gitignore b/ocaml/.gitignore new file mode 100644 index 00000000..059e1de2 --- /dev/null +++ b/ocaml/.gitignore @@ -0,0 +1,5 @@ +ezfio.ml +*.byte +*.native +_build + diff --git a/ocaml/Makefile b/ocaml/Makefile index cbc4d1a3..b2783eb9 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -1,19 +1,37 @@ +# Check if QPACKAGE_ROOT is defined + +ifndef QPACKAGE_ROOT +$(info -------------------- Error --------------------) +$(info QPACKAGE_ROOT undefined. Source the quantum_package.rc script) +$(info -----------------------------------------------) +$(error ) +endif + LIBS= PKGS= OCAMLCFLAGS=-g OCAMLBUILD=ocamlbuild -cflags $(OCAMLCFLAGS) -lflags -g +MLFILES=$(wildcard *.ml) ezfio.ml +MLIFILES=$(wildcard *.mli) +ALL_TESTS=$(patsubst %.ml,%.byte,$(wildcard test_*.ml)) -test_elements.byte: +default: $(ALL_TESTS) -%.inferred.mli: $(wildcard *.ml) +%.inferred.mli: $(MLFILES) $(OCAMLBUILD) $*.inferred.mli -cflags -i -use-ocamlfind $(PKGS) -%.byte: $(wildcard *.ml) $(wildcard *.mli) +%.byte: $(MLFILES) $(MLIFILES) $(OCAMLBUILD) $*.byte -use-ocamlfind $(PKGS) -%.native: $(wildcard *.ml) $(wildcard *.mli) +%.native: $(MLFILES) $(MLIFILES) $(OCAMLBUILD) $*.native -use-ocamlfind $(PKGS) +ezfio.ml: ${QPACKAGE_ROOT}/EZFIO/Ocaml/ezfio.ml + cp ${QPACKAGE_ROOT}/EZFIO/Ocaml/ezfio.ml . + +${QPACKAGE_ROOT}/EZFIO/Ocaml/ezfio.ml: + $(MAKE) -C ${QPACKAGE_ROOT}/src ezfio + clean: rm -rf _build diff --git a/ocaml/README.rst b/ocaml/README.rst new file mode 100644 index 00000000..414428dd --- /dev/null +++ b/ocaml/README.rst @@ -0,0 +1,6 @@ +=============== +Ocaml scripts +=============== + +This directory contains all the scripts that control the input/output +with the user. diff --git a/ocaml/atom.ml b/ocaml/atom.ml index b21b6a7d..1fd9e008 100644 --- a/ocaml/atom.ml +++ b/ocaml/atom.ml @@ -35,6 +35,12 @@ let of_string s = charge = Charge.of_string charge ; coord = Point3d.of_string (String.concat [x; y; z] ?sep:(Some " ")) } + | [ 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:(Some " ")) + } | _ -> raise (AtomError s) ;; @@ -42,6 +48,6 @@ let to_string a = [ Element.to_string a.element ; Charge.to_string a.charge ; Point3d.to_string a.coord ] - |> String.concat ?sep:(Some " ") + |> String.concat ?sep:(Some " ") ;; diff --git a/ocaml/molecule.ml b/ocaml/molecule.ml new file mode 100644 index 00000000..b4845f2c --- /dev/null +++ b/ocaml/molecule.ml @@ -0,0 +1,101 @@ +open Core.Std ;; +open Qptypes ;; + +exception MultiplicityError of string;; + +type t = { + nuclei : Atom.t list ; + elec_alpha : Positive_int.t ; + elec_beta : Positive_int.t ; +} + +let 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 + | [] -> 0. + in + nucl_charge nuclei -. (Float.of_int result) +;; + +let multiplicity m = + Multiplicity.of_alpha_beta m.elec_alpha m.elec_beta +;; + +let name m = + let cm = Float.to_int (charge m) in + let c = + match cm with + | 0 -> "" + | 1 -> " (+)" + | (-1) -> " (-)" + | i when i>1 -> Printf.sprintf " (%d+)" i + | i -> Printf.sprintf " (%d-)" (-i) + in + let mult = Multiplicity.to_string (multiplicity m) in + let { nuclei ; elec_alpha ; elec_beta } = m in + let rec build_list accu = function + | a::rest -> + begin + let e = a.Atom.element in + match (List.Assoc.find accu e) with + | None -> build_list (List.Assoc.add accu e 1) rest + | Some i -> build_list (List.Assoc.add accu e (i+1)) rest + end + | [] -> accu + in + let rec build_name accu = function + | (a, n)::rest -> + let a = Element.to_string a in + begin + match n with + | 1 -> build_name (a::accu) rest + | i when i>1 -> + let tmp = Printf.sprintf "%s%d" a i + in build_name (tmp::accu) rest + | _ -> assert false + end + | [] -> accu + in + let result = build_list [] nuclei |> build_name [c ; ", " ; mult] + in + String.concat (result) +;; + +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) + |> String.concat ?sep:(Some "\n") +;; + +let of_xyz_string s charge' multiplicity' = + let l = String.split s ~on:'\n' + |> List.filter ~f:(fun x -> x <> "") + |> List.map ~f:Atom.of_string + in + let ne = ( charge { + nuclei=l ; + elec_alpha=(Positive_int.of_int 0) ; + elec_beta=(Positive_int.of_int 0) } + |> Float.to_int + )- charge' + |> Positive_int.of_int + in + let (na,nb) = Multiplicity.to_alpha_beta ne multiplicity' in + let result = + { nuclei = l ; + elec_alpha = (Positive_int.of_int na) ; + elec_beta = (Positive_int.of_int nb) } + in + if ((multiplicity result) <> multiplicity') then + let msg = Printf.sprintf + "With %d electrons multiplicity %d is impossible" + (Positive_int.to_int ne) + (Multiplicity.to_int multiplicity') + in + raise (MultiplicityError msg); + else () ; + result +;; diff --git a/ocaml/multiplicity.ml b/ocaml/multiplicity.ml new file mode 100644 index 00000000..89246dec --- /dev/null +++ b/ocaml/multiplicity.ml @@ -0,0 +1,35 @@ +open Qptypes ;; + +type t = Strictly_positive_int.t;; +let of_int = Strictly_positive_int.of_int ;; +let to_int = Strictly_positive_int.to_int ;; + +let to_string m = + match (to_int m) with + | 1 -> "Singlet" + | 2 -> "Doublet" + | 3 -> "Triplet" + | 4 -> "Quartet" + | 5 -> "Quintet" + | 6 -> "Sextet" + | 7 -> "Septet" + | 8 -> "Octet" + | 9 -> "Nonet" + | i -> Printf.sprintf "%d-et" i +;; + +let of_alpha_beta a b = + let a = Positive_int.to_int a + and b = Positive_int.to_int b + in + assert (a >= b); + of_int (1 + a - b) +;; + +let to_alpha_beta ne m = + let ne = Positive_int.to_int ne in + let nb = (ne-(to_int m)+1)/2 in + let na = ne - nb in + assert (na >= nb) ; + (na,nb) +;; diff --git a/ocaml/point3d.ml b/ocaml/point3d.ml index e51003e2..58b39d45 100644 --- a/ocaml/point3d.ml +++ b/ocaml/point3d.ml @@ -30,6 +30,6 @@ let distance p1 p2 = sqrt (distance2 p1 p2) let to_string p = let { x=x ; y=y ; z=z } = p in - Printf.sprintf "%f %f %f" x y z + Printf.sprintf "%f %f %f" x y z ;; diff --git a/ocaml/qptypes.ml b/ocaml/qptypes.ml index da712868..62f00e8f 100644 --- a/ocaml/qptypes.ml +++ b/ocaml/qptypes.ml @@ -5,7 +5,7 @@ module Positive_float : sig end = struct type t = float let to_float x = x - let of_float x = ( assert (x > 0.) ; x ) + let of_float x = ( assert (x >= 0.) ; x ) end @@ -16,7 +16,7 @@ module Strictly_positive_float : sig end = struct type t =float let to_float x = x - let of_float x = ( assert (x >= 0.) ; x ) + let of_float x = ( assert (x > 0.) ; x ) end @@ -27,7 +27,7 @@ module Positive_int : sig end = struct type t = int let to_int x = x - let of_int x = ( assert (x > 0) ; x ) + let of_int x = ( assert (x >= 0) ; x ) end @@ -38,7 +38,7 @@ module Strictly_positive_int : sig end = struct type t = int let to_int x = x - let of_int x = ( assert (x >= 0) ; x ) + let of_int x = ( assert (x > 0) ; x ) end diff --git a/ocaml/test_elements.ml b/ocaml/test_elements.ml index 676e9835..1e91e844 100644 --- a/ocaml/test_elements.ml +++ b/ocaml/test_elements.ml @@ -1,6 +1,6 @@ let test_module () = - let atom = Elements.of_string "Cobalt" in - Printf.printf "%s %d\n" (Elements.to_string atom) (Elements.charge atom) + let atom = Element.of_string "Cobalt" in + Printf.printf "%s %d\n" (Element.to_string atom) (Element.charge atom) ;; test_module ();; diff --git a/ocaml/test_molecule.ml b/ocaml/test_molecule.ml new file mode 100644 index 00000000..d648e9ce --- /dev/null +++ b/ocaml/test_molecule.ml @@ -0,0 +1,35 @@ +open Core.Std ;; +open Qptypes ;; + +let test_molecule () = + let xyz = +" +H 1.0 0.54386314 0.00000000 -3.78645152 +O 8.0 1.65102147 0.00000000 -2.35602344 +H 1.0 0.54386314 0.00000000 -0.92559535 +" + in + + try ignore (Molecule.of_xyz_string xyz 1 (Strictly_positive_int.of_int 1)) + with + | Molecule.MultiplicityError _ -> print_string "OK\n" + ; + print_string "---\n"; + let m = Molecule.of_xyz_string xyz 0 (Strictly_positive_int.of_int 1) + in print_endline (Molecule.name m) ; + let m = Molecule.of_xyz_string xyz 1 (Strictly_positive_int.of_int 2) + in print_endline (Molecule.name m) ; + + let xyz = +" +H 0.54386314 0.00000000 -3.78645152 +O 1.65102147 0.00000000 -2.35602344 +H 0.54386314 0.00000000 -0.92559535 +" + in + let m = Molecule.of_xyz_string xyz (-2) (Strictly_positive_int.of_int 1) + in print_endline (Molecule.name m) ; + print_string (Molecule.to_string m); +;; + +test_molecule ();; diff --git a/src/Makefile.common b/src/Makefile.common index 86c96fbb..0aa47245 100644 --- a/src/Makefile.common +++ b/src/Makefile.common @@ -96,6 +96,7 @@ $(EZFIO): $(wildcard $(QPACKAGE_ROOT)/src/*.ezfio_config) $(wildcard $(QPACKAGE_ @cp $(wildcard $(QPACKAGE_ROOT)/src/*.ezfio_config) $(wildcard $(QPACKAGE_ROOT)/src/*/*.ezfio_config) $(EZFIO_DIR)/config @cd $(EZFIO_DIR) ; export FC="$(FC)" ; export FCFLAGS="$(FCFLAGS)" ; export IRPF90="$(IRPF90)" ; $(MAKE) ; $(MAKE) Python +ezfio: $(EZFIO) # Create symbolic links of other modules diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 0ef6daa9..e69de29b 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1,2 +0,0 @@ -AOs BiInts Bitmask CISD CISD_SC2_selected CISD_selected DensityMatrix Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Perturbation SC2 Selectors_full SingleRefMethod Utils -