From 4696e3d5a9cddf48ea45944f4c6def6e4782b07a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Jan 2016 15:44:15 +0100 Subject: [PATCH] Pseudo in OCaml --- data/pseudo/bfd | 552 ++++++++++++++++++++++++++++++ ocaml/.gitignore | 2 + ocaml/Gto.ml | 12 +- ocaml/Makefile | 13 +- ocaml/Pseudo.ml | 273 +++++++++++++++ ocaml/qp_create_ezfio_from_xyz.ml | 264 ++++++++++++-- ocaml/qp_edit.ml | 24 +- ocaml/qptypes_generator.ml | 6 +- 8 files changed, 1090 insertions(+), 56 deletions(-) create mode 100644 data/pseudo/bfd create mode 100644 ocaml/Pseudo.ml diff --git a/data/pseudo/bfd b/data/pseudo/bfd new file mode 100644 index 00000000..8dfaffaf --- /dev/null +++ b/data/pseudo/bfd @@ -0,0 +1,552 @@ +H GEN 0 0 +3 +1.00000000 1 4.47692410 +4.47692410 3 2.97636451 +-4.32112340 2 3.01841596 + +He GEN 0 0 +3 +2.00000000 1 9.84368811 +19.68737622 3 10.94516233 +-17.20570338 2 12.03715264 + +Li GEN 2 1 +3 +1.00000000 1 5.41040609 +5.41040609 3 2.70520138 +-4.60151975 2 2.07005488 +1 +7.09172172 2 1.34319829 + +Be GEN 2 1 +3 +2.00000000 1 4.58686001 +9.17372003 3 2.29371778 +-8.12599146 2 2.10401964 +1 +14.90499810 2 2.71723988 + +B GEN 2 1 +3 +3.00000000 1 5.40423964 +16.21271892 3 5.71678458 +-11.86640633 2 4.48974455 +1 +15.49737620 2 3.43781634 + +C GEN 2 1 +3 +4.00000000 1 8.35973821 +33.43895285 3 4.48361888 +-19.17537323 2 3.93831258 +1 +22.55164191 2 5.02991637 + +N GEN 2 1 +3 +5.00000000 1 9.23501007 +46.17505034 3 7.66830008 +-30.18893534 2 7.34486070 +1 +31.69720409 2 6.99536540 + +O GEN 2 1 +3 +6.00000000 1 9.29793903 +55.78763416 3 8.86492204 +-38.81978498 2 8.62925665 +1 +38.41914135 2 8.71924452 + +F GEN 2 1 +3 +7.00000000 1 11.39210685 +79.74474797 3 10.74911370 +-49.45159098 2 10.45120693 +1 +50.25646328 2 11.30345826 + +Ne GEN 2 1 +3 +8.00000000 1 10.74945199 +85.99561593 3 10.19801460 +-56.79004456 2 10.18694048 +1 +55.11144535 2 12.85042963 + +Na GEN 10 2 +3 +1.00000000 1 5.35838717 +5.35838717 3 3.67918975 +-2.07764789 2 1.60507673 +1 +10.69640234 2 1.32389367 +1 +10.11238853 2 1.14052020 + +Mg GEN 10 2 +3 +2.00000000 1 4.48537898 +8.97075796 3 2.24268949 +-7.72153408 2 1.59710474 +1 +15.07848048 2 1.57188656 +1 +12.37888383 2 1.42757357 + +Al GEN 10 2 +3 +3.00000000 1 3.07301275 +9.21903825 3 9.97055633 +-9.65888637 2 2.64134660 +1 +17.16680112 2 1.87284747 +1 +14.22120694 2 1.79397208 + +Si GEN 10 2 +3 +4.00000000 1 1.80721061 +7.22884246 3 9.99633089 +-13.06725590 2 2.50043232 +1 +21.20531613 2 2.26686403 +1 +15.43693603 2 2.11659661 + +P GEN 10 2 +3 +5.00000000 1 2.02622810 +10.13114051 3 9.95970113 +-14.94375088 2 2.74841795 +1 +23.62479480 2 2.60470698 +1 +18.18547203 2 2.54957900 + +S GEN 10 2 +3 +6.00000000 1 2.42178462 +14.53070769 3 6.74148698 +-17.52965289 2 3.06094751 +1 +25.99260928 2 2.94272173 +1 +18.93356489 2 2.84566981 + +Cl GEN 10 2 +3 +7.00000000 1 2.41079533 +16.87556731 3 5.29139158 +-18.80917558 2 2.91105513 +1 +28.59107316 2 3.34528827 +1 +19.37583724 2 3.12408551 + +Ar GEN 10 2 +3 +8.00000000 1 3.09403094 +24.75224749 3 6.53700323 +-20.38446872 2 3.35769859 +1 +30.67006675 2 3.68203169 +1 +20.84338017 2 3.45735664 + +K GEN 10 2 +3 +9.00000000 1 1.61692238 +14.55230143 3 2.94781323 +-28.90972088 2 2.34518380 +1 +85.63621269 2 6.20715645 +1 +30.81028404 2 4.43761067 + +Ca GEN 10 2 +3 +10.00000000 1 1.83448200 +18.34481997 3 2.93139337 +-33.15102838 2 2.50808985 +1 +90.63863403 2 6.91765747 +1 +31.99145044 2 4.89613174 + +Sc GEN 10 2 +3 +11.00000000 1 1.64916555 +18.14082106 3 3.06833288 +-35.19310566 2 2.41985389 +1 +97.62913482 2 7.60319353 +1 +33.97033635 2 5.31121835 + +Ti GEN 10 2 +3 +12.00000000 1 1.85755224 +22.29062683 3 3.30638246 +-41.26280518 2 2.70879079 +1 +96.94231089 2 8.03040157 +1 +38.01641313 2 5.93291405 + +V GEN 10 2 +3 +13.00000000 1 2.16361754 +28.12702797 3 4.07901780 +-48.27656329 2 3.21436396 +1 +96.23226580 2 8.44326050 +1 +41.58043539 2 6.53136059 + +Cr GEN 10 2 +3 +14.00000000 1 2.94337662 +41.20727267 3 3.40750349 +-55.51413840 2 3.33587110 +1 +103.26274052 2 9.14231779 +1 +45.80124572 2 7.21220771 + +Mn GEN 10 2 +3 +15.00000000 1 3.29524605 +49.42869068 3 3.61599152 +-61.66925169 2 3.57405761 +1 +112.85037953 2 9.99154195 +1 +49.18832867 2 7.80925218 + +Fe GEN 10 2 +3 +16.00000000 1 3.72075632 +59.53210107 3 3.92321272 +-68.75847841 2 3.89595440 +1 +112.92561163 2 10.42343546 +1 +52.55882759 2 8.41664076 + +Co GEN 10 2 +3 +17.00000000 1 4.18819322 +71.19928469 3 4.42968808 +-77.65278252 2 4.39800669 +1 +113.90484511 2 10.86075441 +1 +56.10698766 2 9.05202771 + +Ni GEN 10 2 +3 +18.00000000 1 4.22878924 +76.11820629 3 4.46202418 +-82.85330412 2 4.44960456 +1 +120.84003628 2 11.62700064 +1 +58.54148726 2 9.57327085 + +Cu GEN 10 2 +3 +19.00000000 1 6.25149628 +118.77842940 3 6.72725326 +-105.89982403 2 6.61024592 +1 +127.35069424 2 12.36568715 +1 +71.22984900 2 11.16072939 + +Zn GEN 10 2 +3 +20.00000000 1 5.25282726 +105.05654526 3 5.85433298 +-105.08806248 2 5.80452897 +1 +123.87006927 2 12.52174964 +1 +72.33499364 2 11.56019052 + +Ga GEN 28 3 +3 +3.00000000 1 1.22544253 +3.67632759 3 5.71065133 +-11.23828733 2 1.33931968 +1 +57.88512249 2 2.48772664 +1 +43.67871044 2 2.12489462 +1 +17.97137628 2 1.27124173 + +Ge GEN 28 3 +3 +4.00000000 1 1.88492329 +7.53969315 3 9.98137268 +-13.13622589 2 1.98008364 +1 +61.26369269 2 3.03315885 +1 +55.52495744 2 2.76564031 +1 +23.49168485 2 1.66026543 + +As GEN 28 3 +3 +5.00000000 1 1.21796059 +6.08980295 3 9.96302171 +-14.92625816 2 1.86158567 +1 +73.75553709 2 3.54546456 +1 +68.03580909 2 3.28664249 +1 +23.32540737 2 1.95862616 + +Se GEN 28 3 +3 +6.00000000 1 1.73494732 +10.40968393 3 6.91632737 +-17.83199463 2 3.10551681 +1 +85.94238004 2 4.67503354 +1 +78.84838432 2 4.34256579 +1 +30.92151589 2 2.61905005 + +Br GEN 28 3 +3 +7.00000000 1 1.86793881 +13.07557164 3 5.30536536 +-18.79056037 2 3.32134623 +1 +88.58537968 2 5.17694821 +1 +79.43718432 2 4.80714881 +1 +29.35463757 2 3.03534088 + +Kr GEN 28 3 +3 +8.00000000 1 1.72397711 +13.79181690 3 6.70510242 +-22.77215308 2 2.75463922 +1 +92.78570269 2 4.85045356 +1 +80.37767796 2 4.52350391 +1 +31.36172413 2 3.04556109 + +Rb GEN 28 3 +3 +9.00000000 1 1.93326235 +17.39936111 3 3.14473834 +-26.53731001 2 2.67367304 +1 +85.71767228 2 5.15576814 +1 +48.02236719 2 4.13480924 +1 +32.98704385 2 3.52393676 + +Sr GEN 28 3 +3 +10.00000000 1 1.51098327 +15.10983267 3 4.08312080 +-26.67485196 2 2.20539101 +1 +91.10888571 2 4.83388773 +1 +51.70644278 2 3.81526734 +1 +31.98712514 2 3.20797541 + +In GEN 46 3 +3 +3.00000000 1 0.84508372 +2.53525117 3 5.66019931 +-7.66579852 2 0.85764327 +1 +58.16918845 2 1.87837596 +1 +43.63891951 2 1.51547534 +1 +17.93363847 2 0.83917399 + +Sn GEN 46 3 +3 +4.00000000 1 0.61334103 +2.45336413 3 5.67826592 +-9.33070594 2 0.86945138 +1 +58.04190484 2 2.01380769 +1 +43.68447157 2 1.63883815 +1 +17.95660523 2 0.92346533 + +Sb GEN 46 3 +3 +5.00000000 1 0.74773404 +3.73867018 3 5.79307847 +-13.88202267 2 1.07909250 +1 +57.56076138 2 2.04356327 +1 +43.88817098 2 1.70062095 +1 +17.82275877 2 1.00414410 + +Te GEN 46 3 +3 +6.00000000 1 1.28872858 +7.73237146 3 5.81046154 +-13.65642757 2 1.55436985 +1 +57.52907325 2 2.39157624 +1 +43.70165092 2 1.96781367 +1 +18.64325026 2 1.13849722 + +I GEN 46 3 +3 +7.00000000 1 1.01780923 +7.12466460 3 3.60136147 +-29.18419372 2 1.68345378 +1 +108.68417388 2 2.80278521 +1 +99.35380694 2 2.51494051 +1 +41.32653157 2 1.56672279 + +Xe GEN 46 3 +3 +8.00000000 1 1.22029027 +9.76232213 3 5.14625292 +-32.16318995 2 2.10563577 +1 +110.93633943 2 3.20320603 +1 +99.49007680 2 2.86062806 +1 +41.81892440 2 1.74523568 + +Cs GEN 46 3 +3 +9.00000000 1 1.42083294 +12.78749648 3 3.04497432 +-21.48605417 2 3.27432859 +1 +87.71882150 2 4.74273103 +1 +50.11684516 2 3.50080269 +1 +34.28999477 2 2.05168748 + +Ba GEN 46 3 +3 +10.00000000 1 1.75114866 +17.51148658 3 7.10519292 +-27.30238921 2 1.92434265 +1 +86.12190858 2 2.79769987 +1 +52.48439280 2 2.16940023 +1 +41.96330733 2 1.76672762 + +Tl GEN 78 4 +3 +3.00000000 1 1.34540299 +4.03620897 3 1.80622564 +-21.68751419 2 0.88272932 +1 +91.69120822 2 1.55008480 +1 +85.52036047 2 1.34656124 +1 +62.00124014 2 0.92193030 +1 +45.59483689 2 0.95782546 + +Pb GEN 78 4 +3 +4.00000000 1 1.33905502 +5.35622008 3 3.57680327 +-25.11165802 2 1.08584447 +1 +121.69447681 2 1.89957262 +1 +114.36466627 2 1.64009233 +1 +49.32959048 2 0.93051806 +1 +45.59434323 2 1.07638351 + +Bi GEN 78 4 +3 +5.00000000 1 1.32337887 +6.61689433 3 3.57500194 +-25.11872438 2 1.26988106 +1 +121.69346942 2 2.17897203 +1 +114.36404033 2 1.85248885 +1 +59.32816833 2 1.10815262 +1 +55.59430794 2 1.13749297 + +Po GEN 78 4 +3 +6.00000000 1 1.12915012 +6.77490072 3 3.57645639 +-27.15687722 2 1.43566151 +1 +123.68767392 2 2.45166968 +1 +116.36047985 2 2.06261668 +1 +81.32822758 2 1.28212112 +1 +80.59441478 2 1.32883356 + +At GEN 78 4 +3 +7.00000000 1 1.03089969 +7.21629781 3 3.69262441 +-29.19452905 2 1.45693225 +1 +128.68209352 2 2.56609159 +1 +119.35613288 2 2.15280178 +1 +80.33152157 2 1.37630746 +1 +80.58962589 2 1.47835294 + +Rn GEN 78 4 +3 +8.00000000 1 1.09138091 +8.73104728 3 3.03187023 +-29.28994938 2 1.65026568 +1 +128.63606138 2 2.92260183 +1 +119.39674828 2 2.42144059 +1 +80.33096401 2 1.49781359 +1 +80.58961959 2 1.47976575 + diff --git a/ocaml/.gitignore b/ocaml/.gitignore index 2ad52540..a655960d 100644 --- a/ocaml/.gitignore +++ b/ocaml/.gitignore @@ -50,6 +50,8 @@ test_molecule test_molecule.byte test_point3d test_point3d.byte +test_pseudo +test_pseudo.byte test_queuing_system test_queuing_system.byte test_task_server diff --git a/ocaml/Gto.ml b/ocaml/Gto.ml index e3c5d62a..69aeba37 100644 --- a/ocaml/Gto.ml +++ b/ocaml/Gto.ml @@ -1,5 +1,5 @@ -open Core.Std;; -open Qptypes;; +open Core.Std +open Qptypes exception GTO_Read_Failure of string exception End_Of_Basis @@ -8,7 +8,7 @@ type t = { sym : Symmetry.t ; lc : ((Primitive.t * AO_coef.t) list) } with sexp -;; + let of_prim_coef_list pc = let (p,c) = List.hd_exn pc in @@ -27,7 +27,7 @@ let of_prim_coef_list pc = { sym = sym ; lc = pc } -;; + let read_one in_channel = @@ -65,7 +65,7 @@ let read_one in_channel = in read_lines [] n |> List.rev |> of_prim_coef_list -;; + (** Transform the gto to a string *) @@ -86,5 +86,5 @@ let to_string { sym = sym ; lc = lc } = in (do_work [result] 1 lc) |> String.concat ~sep:"\n" -;; + diff --git a/ocaml/Makefile b/ocaml/Makefile index 01473ec2..31330c66 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -13,17 +13,18 @@ LIBS= PKGS= OCAMLCFLAGS="-g -warn-error A" OCAMLBUILD=ocamlbuild -j 0 -syntax camlp4o -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) -MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml Input_auto_generated.ml qp_edit.ml -MLIFILES=$(wildcard *.mli) +MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml Input_auto_generated.ml qp_edit.ml +MLIFILES=$(wildcard *.mli) git ALL_TESTS=$(patsubst %.ml,%.byte,$(wildcard test_*.ml)) ALL_EXE=$(patsubst %.ml,%.native,$(wildcard qp_*.ml)) qp_edit.native .PHONY: executables default remake_executables -default: $(ALL_TESTS) $(ALL_EXE) .gitignore +default: $(ALL_EXE) .gitignore +tests: $(ALL_TESTS) -.gitignore: $(MLFILES) +.gitignore: $(MLFILES) $(MLIFILES) @for i in .gitignore ezfio.ml Qptypes.ml Git.ml qptypes_generator.byte _build $(ALL_EXE) $(ALL_TESTS) \ $(patsubst %.ml,%,$(wildcard test_*.ml)) $(patsubst %.ml,%,$(wildcard qp_*.ml)) \ $(shell grep Input Input_auto_generated.ml | awk '{print $$2 ".ml"}') \ @@ -69,8 +70,10 @@ ezfio.ml: ${QP_ROOT}/install/EZFIO/Ocaml/ezfio.ml qptypes_generator.byte: qptypes_generator.ml $(OCAMLBUILD) qptypes_generator.byte -use-ocamlfind -Git.ml Qptypes.ml: qptypes_generator.byte +Qptypes.ml: qptypes_generator.byte ./qptypes_generator.byte > Qptypes.ml + +git: ./create_git_sha1.sh ${QP_ROOT}/install/EZFIO/Ocaml/ezfio.ml: diff --git a/ocaml/Pseudo.ml b/ocaml/Pseudo.ml new file mode 100644 index 00000000..3fb4736e --- /dev/null +++ b/ocaml/Pseudo.ml @@ -0,0 +1,273 @@ +open Core.Std +open Qptypes + + +module Primitive_local : sig + + type t = { + expo : AO_expo.t ; + r_power : R_power.t ; + } with sexp + + val of_expo_r_power : AO_expo.t -> R_power.t -> t + val to_string : t -> string + +end = struct + + type t = { + expo : AO_expo.t ; + r_power : R_power.t ; + } with sexp + + let of_expo_r_power dz n = + { expo = dz ; r_power = n } + + let to_string p = + Printf.sprintf "(%d, %f)" + (R_power.to_int p.r_power) + (AO_expo.to_float p.expo) +end + + +module Primitive_non_local : sig + + type t = { + expo : AO_expo.t ; + r_power : R_power.t ; + proj : Positive_int.t + } with sexp + + val of_proj_expo_r_power : Positive_int.t -> AO_expo.t -> R_power.t -> t + val to_string : t -> string + +end = struct + + type t = { + expo : AO_expo.t ; + r_power : R_power.t ; + proj : Positive_int.t + } with sexp + + let of_proj_expo_r_power p dz n = + { expo = dz ; r_power = n ; proj = p } + + let to_string p = + Printf.sprintf "(%d, %f, %d)" + (R_power.to_int p.r_power) + (AO_expo.to_float p.expo) + (Positive_int.to_int p.proj) +end + + + + +type t = { + element : Element.t ; + n_elec : Positive_int.t ; + local : (Primitive_local.t * AO_coef.t ) list ; + non_local : (Primitive_non_local.t * AO_coef.t ) list +} with sexp + +let empty e = + { element = e; + n_elec = Positive_int.of_int 0; + local = []; + non_local = []; + } + +(** Transform the local component of the pseudopotential to a string *) +let to_string_local = function +| [] -> "" +| t -> + "Local component:" :: + ( Printf.sprintf "%20s %8s %20s" "Coeff." "r^n" "Exp." ) :: + ( List.map t ~f:(fun (l,c) -> Printf.sprintf "%20f %8d %20f" + (AO_coef.to_float c) + (R_power.to_int l.Primitive_local.r_power) + (AO_expo.to_float l.Primitive_local.expo) + ) ) + |> String.concat ~sep:"\n" + + +(** Transform the non-local component of the pseudopotential to a string *) +let to_string_non_local = function +| [] -> "" +| t -> + "Non-local component:" :: + ( Printf.sprintf "%20s %8s %20s %8s" "Coeff." "r^n" "Exp." "Proj") :: + ( List.map t ~f:(fun (l,c) -> + let p = + Positive_int.to_int l.Primitive_non_local.proj + in + Printf.sprintf "%20f %8d %20f |%d><%d|" + (AO_coef.to_float c) + (R_power.to_int l.Primitive_non_local.r_power) + (AO_expo.to_float l.Primitive_non_local.expo) + p p + ) ) + |> String.concat ~sep:"\n" + +(** Transform the Pseudopotential to a string *) +let to_string t = + + Printf.sprintf "%s %d electrons removed" + (Element.to_string t.element) + (Positive_int.to_int t.n_elec) + :: to_string_local t.local + :: to_string_non_local t.non_local + :: [] + |> List.filter ~f:(fun x -> x <> "") + |> String.concat ~sep:"\n" + + +(** Find an element in the file *) +let find in_channel element = + In_channel.seek in_channel 0L; + + let element_read, old_pos = + ref Element.X, + ref (In_channel.pos in_channel) + in + while !element_read <> element + do + let buffer = + old_pos := In_channel.pos in_channel; + match In_channel.input_line in_channel with + | Some line -> String.split ~on:' ' line + |> List.hd_exn + | None -> "" + in + try + element_read := Element.of_string buffer + with + | Element.ElementError _ -> () + done ; + In_channel.seek in_channel !old_pos; + !element_read + + +(** Read the Pseudopotential in GAMESS format *) +let read_element in_channel element = + ignore (find in_channel element); + + let rec read result = + match In_channel.input_line in_channel with + | None -> result + | Some line -> + if (String.strip line = "") then + result + else + read (line::result) + in + + let data = + read [] + |> List.rev + in + + let debug_data = + String.concat ~sep:"\n" data + in + + let decode_first_line = function + | first_line :: rest -> + begin + let first_line_split = + String.split first_line ~on:' ' + |> List.filter ~f:(fun x -> (String.strip x) <> "") + in + match first_line_split with + | e :: "GEN" :: n :: p -> + { element = Element.of_string e ; + n_elec = Int.of_string n |> Positive_int.of_int ; + local = [] ; + non_local = [] + }, rest + | _ -> failwith ( + Printf.sprintf "Unable to read Pseudopotential : \n%s\n" + debug_data ) + end + | _ -> failwith ("Error reading pseudopotential\n"^debug_data) + in + + let rec loop create_primitive accu = function + | (0,rest) -> List.rev accu, rest + | (n,line::rest) -> + begin + match + String.split line ~on:' ' + |> List.filter ~f:(fun x -> String.strip x <> "") + with + | c :: i :: e :: [] -> + let i = + Int.of_string i + in + let elem = + ( create_primitive + (Float.of_string e |> AO_expo.of_float) + (i-2 |> R_power.of_int), + Float.of_string c |> AO_coef.of_float + ) + in + loop create_primitive (elem::accu) (n-1, rest) + | _ -> failwith ("Error reading pseudopotential\n"^debug_data) + end + | _ -> failwith ("Error reading pseudopotential\n"^debug_data) + in + + let decode_local (pseudo,data) = + let decode_local_n n rest = + let result, rest = + loop Primitive_local.of_expo_r_power [] (Positive_int.to_int n,rest) + in + { pseudo with local = result }, rest + in + match data with + | n :: rest -> + let n = + String.strip n + |> Int.of_string + |> Positive_int.of_int + in + decode_local_n n rest + | _ -> failwith ("Unable to read (non-)local pseudopotential\n"^debug_data) + in + + let decode_non_local (pseudo,data) = + let decode_non_local_n proj n (pseudo,data) = + let result, rest = + loop (Primitive_non_local.of_proj_expo_r_power proj) + [] (Positive_int.to_int n, data) + in + { pseudo with non_local = pseudo.non_local @ result }, rest + in + let rec new_proj (pseudo,data) proj = + match data with + | n :: rest -> + let n = + String.strip n + |> Int.of_string + |> Positive_int.of_int + in + let result = + decode_non_local_n proj n (pseudo,rest) + and proj_next = + (Positive_int.to_int proj)+1 + |> Positive_int.of_int + in + new_proj result proj_next + | _ -> pseudo + in + new_proj (pseudo,data) (Positive_int.of_int 0) + in + + decode_first_line data + |> decode_local + |> decode_non_local + + + + +include To_md5 +let to_md5 = to_md5 sexp_of_t + diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index ee7e7d40..3a1dda1a 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -15,11 +15,12 @@ let spec = ~doc:"float Add dummy atoms. x * (covalent radii of the atoms)" +> flag "m" (optional_with_default 1 int) ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1." - +> flag "p" no_arg - ~doc:" Using pseudopotentials" + +> flag "p" (optional string) + ~doc:"string Name of the pseudopotential" +> anon ("xyz_file" %: file ) +(** Handle dummy atoms placed on bonds *) let dummy_centers ~threshold ~molecule ~nuclei = let d = Molecule.distance_matrix molecule @@ -68,6 +69,7 @@ let dummy_centers ~threshold ~molecule ~nuclei = ) +(** Returns the list of available basis sets *) let list_basis () = let basis_list = Qpackage.root ^ "/install/emsl/EMSL_api.py list_basis" @@ -84,6 +86,7 @@ let list_basis () = |> String.concat ~sep:"\t" +(** Run the program *) let run ?o b c d m p xyz_file = (* Read molecule *) @@ -94,16 +97,15 @@ let run ?o b c d m p xyz_file = let dummy = dummy_centers ~threshold:d ~molecule ~nuclei:molecule.Molecule.nuclei in -(* - List.iter dummy ~f:(fun x -> - Printf.printf "%s\n" (Atom.to_string ~units:Units.Angstrom x) - ); -*) let nuclei = molecule.Molecule.nuclei @ dummy in + (********** + Basis set + **********) + let basis_table = Hashtbl.Poly.create () in @@ -133,11 +135,7 @@ let run ?o b c d m p xyz_file = let fetch_channel basis = let command = - if (p) then - Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename - ^ "." ^ basis ^ "\" \"" ^ basis ^"\" pseudo" - else - Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename + Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename ^ "." ^ basis ^ "\" \"" ^ basis ^"\"" in match Sys.is_file basis with @@ -206,6 +204,89 @@ let run ?o b c d m p xyz_file = |> List.rev_map ~f:String.strip |> build_basis; + + + (*************** + Pseudopotential + ***************) + + let pseudo_table = + Hashtbl.Poly.create () + in + + (* Open pseudo channels *) + let pseudo_channel element = + let key = + Element.to_string element + in + Hashtbl.find pseudo_table key + in + let temp_filename = + Filename.temp_file "qp_create_" ".pseudo" + in + let () = + Sys.remove temp_filename + in + + let fetch_channel pseudo = + match Sys.is_file pseudo with + | `Yes -> In_channel.create pseudo + | _ -> failwith ("Pseudo file "^pseudo^" not found.") + in + + let rec build_pseudo = function + | [] -> () + | elem_and_pseudo_name :: rest -> + begin + match (String.lsplit2 ~on:':' elem_and_pseudo_name) with + | None -> (* Principal pseudo *) + begin + let pseudo = + elem_and_pseudo_name + in + let new_channel = + fetch_channel pseudo + in + List.iter nuclei ~f:(fun elem-> + let key = + Element.to_string elem.Atom.element + in + match Hashtbl.add pseudo_table ~key:key ~data:new_channel with + | `Ok -> () + | `Duplicate -> () + ) + end + | Some (key, pseudo) -> (*Aux pseudo *) + begin + let elem = + Element.of_string key + and pseudo = + String.lowercase pseudo + in + let key = + Element.to_string elem + in + let new_channel = + fetch_channel pseudo + in + begin + match Hashtbl.add pseudo_table ~key:key ~data:new_channel with + | `Ok -> () + | `Duplicate -> failwith ("Duplicate definition of pseudo for "^(Element.to_long_string elem)) + end + end + end; + build_pseudo rest + in + let () = + match p with + | None -> () + | Some p -> + String.split ~on:'|' p + |> List.rev_map ~f:String.strip + |> build_pseudo + in + (* Build EZFIO File name *) let ezfio_file = match o with @@ -223,6 +304,44 @@ let run ?o b c d m p xyz_file = (* Create EZFIO *) Ezfio.set_file ezfio_file; + (* Write Pseudo *) + let pseudo = + List.map nuclei ~f:(fun x -> + match pseudo_channel x.Atom.element with + | Some channel -> Pseudo.read_element channel x.Atom.element + | None -> Pseudo.empty x.Atom.element + ) + in + List.iter pseudo ~f:(fun x -> Pseudo.to_string x |> print_endline) ; + + let molecule = + let n_elec_to_remove = + List.fold pseudo ~init:0 ~f:(fun accu x -> + accu + (Positive_int.to_int x.Pseudo.n_elec)) + in + { Molecule.elec_alpha = + (Elec_alpha_number.to_int molecule.Molecule.elec_alpha) + - n_elec_to_remove/2 + |> Elec_alpha_number.of_int; + Molecule.elec_beta = + (Elec_beta_number.to_int molecule.Molecule.elec_beta) + - (n_elec_to_remove - n_elec_to_remove/2) + |> Elec_beta_number.of_int; + Molecule.nuclei = + let charges = + List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec + |> Float.of_int) + |> Array.of_list + in + List.mapi molecule.Molecule.nuclei ~f:(fun i x -> + { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i) + |> Charge.of_float } + ) + } + in + List.iter molecule.Molecule.nuclei ~f:(fun x -> print_endline (Atom.to_string x ~units:Units.Bohr)); + + (* Write Electrons *) Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int molecule.Molecule.elec_alpha ) ; @@ -247,6 +366,92 @@ let run ?o b c d m p xyz_file = Ezfio.set_nuclei_nucl_coord (Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords); + + (* Write pseudopotential *) + let () = + match p with + | None -> Ezfio.set_pseudo_do_pseudo false + | _ -> Ezfio.set_pseudo_do_pseudo true + in + + let klocmax = + List.fold pseudo ~init:0 ~f:(fun accu x -> + let x = + List.length x.Pseudo.local + in + if (x > accu) then x + else accu + ) + and kmax = + List.fold pseudo ~init:0 ~f:(fun accu x -> + let x = + List.length x.Pseudo.non_local + in + if (x > accu) then x + else accu + ) + and lmax = + List.fold pseudo ~init:0 ~f:(fun accu x -> + let x = + List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) -> + let x = + Positive_int.to_int x.Pseudo.Primitive_non_local.proj + in + if (x > accu) then x + else accu + ) + in + if (x > accu) then x + else accu + ) + in + + let () = + Ezfio.set_pseudo_pseudo_klocmax klocmax; + Ezfio.set_pseudo_pseudo_kmax kmax; + Ezfio.set_pseudo_pseudo_lmax lmax; + let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k = + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. , + Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0 + in + List.iteri pseudo ~f:(fun j x -> + List.iteri x.Pseudo.local ~f:(fun i (y,c) -> + tmp_array_v_k.(i).(j) <- AO_coef.to_float c; + let y, z = + AO_expo.to_float y.Pseudo.Primitive_local.expo, + R_power.to_int y.Pseudo.Primitive_local.r_power + in + tmp_array_dz_k.(i).(j) <- y; + tmp_array_n_k.(i).(j) <- z; + ) + ); + let v_k = + Array.map tmp_array_v_k ~f:Array.to_list + |> Array.to_list + |> List.concat + in + Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data: v_k + |> Ezfio.set_pseudo_pseudo_v_k ; + let dz_k = + Array.map tmp_array_dz_k ~f:Array.to_list + |> Array.to_list + |> List.concat + in + Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data: dz_k + |> Ezfio.set_pseudo_pseudo_dz_k; + let n_k = + Array.map tmp_array_n_k ~f:Array.to_list + |> Array.to_list + |> List.concat + in + Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data: n_k + |> Ezfio.set_pseudo_pseudo_n_k; + in + + + + (* Write Basis set *) let basis = @@ -296,7 +501,9 @@ let run ?o b c d m p xyz_file = if x > s then x else s) ao_prim_num in - let gtos = List.map long_basis ~f:(fun (_,x,_) -> x) in + let gtos = + List.map long_basis ~f:(fun (_,x,_) -> x) + in let create_expo_coef ec = let coefs = @@ -329,25 +536,18 @@ let run ?o b c d m p xyz_file = 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) ; - 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) ; - - - (* Doesn't work... *) - if (p) then - begin - Qpackage.root ^ "/scripts/pseudo/put_pseudo_in_ezfio.py " ^ ezfio_file - |> Sys.command_exn - end; - + let () = + 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) ; + 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) ; + in match Input.Ao_basis.read () with | None -> failwith "Error in basis" | Some x -> Input.Ao_basis.write x diff --git a/ocaml/qp_edit.ml b/ocaml/qp_edit.ml index c927d3e9..a693aa2f 100644 --- a/ocaml/qp_edit.ml +++ b/ocaml/qp_edit.ml @@ -13,31 +13,31 @@ This file is autogenerad by (** Keywords used to define input sections *) type keyword = | Ao_basis -| Determinants | Determinants_by_hand | Electrons -| Hartree_fock -| Integrals_bielec | Mo_basis | Nuclei +| Determinants +| Integrals_bielec +| Pseudo | Perturbation | Properties -| Pseudo +| Hartree_fock ;; let keyword_to_string = function | Ao_basis -> "AO basis" | Determinants_by_hand -> "Determinants_by_hand" -| Determinants -> "Determinants" | Electrons -> "Electrons" -| Hartree_fock -> "Hartree_fock" -| Integrals_bielec -> "Integrals_bielec" | Mo_basis -> "MO basis" | Nuclei -> "Molecule" +| Determinants -> "Determinants" +| Integrals_bielec -> "Integrals_bielec" +| Pseudo -> "Pseudo" | Perturbation -> "Perturbation" | Properties -> "Properties" -| Pseudo -> "Pseudo" +| Hartree_fock -> "Hartree_fock" ;; @@ -94,10 +94,10 @@ let get s = f Pseudo.(read, to_rst) | Perturbation -> f Perturbation.(read, to_rst) - | Hartree_fock -> - f Hartree_fock.(read, to_rst) | Properties -> f Properties.(read, to_rst) + | Hartree_fock -> + f Hartree_fock.(read, to_rst) end with | Sys_error msg -> (Printf.eprintf "Info: %s\n%!" msg ; "") @@ -139,8 +139,8 @@ let set str s = | Integrals_bielec -> write Integrals_bielec.(of_rst, write) s | Pseudo -> write Pseudo.(of_rst, write) s | Perturbation -> write Perturbation.(of_rst, write) s - | Hartree_fock -> write Hartree_fock.(of_rst, write) s | Properties -> write Properties.(of_rst, write) s + | Hartree_fock -> write Hartree_fock.(of_rst, write) s | Electrons -> write Electrons.(of_rst, write) s | Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s | Nuclei -> write Nuclei.(of_rst, write) s @@ -192,8 +192,8 @@ let run check_only ezfio_filename = Integrals_bielec ; Pseudo ; Perturbation ; - Hartree_fock ; Properties ; + Hartree_fock ; Mo_basis; Determinants_by_hand ; ] diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index 44901af8..f3a5513e 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -79,6 +79,10 @@ let input_data = " * AO_prim_number : int assert (x > 0) ; +* R_power : int + assert (x >= -2) ; + assert (x <= 8) ; + * Threshold : float assert (x >= 0.) ; assert (x <= 1.) ; @@ -138,7 +142,7 @@ let input_ezfio = " * Det_number : int determinants_n_det 1 : 100000000 - More than 100 million determinants + More than 100 million of determinants " ;;