LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-01-03 10:05:57 +01:00
Code Issues Releases Wiki Activity

Pseudo in OCaml

Browse Source
This commit is contained in:
Anthony Scemama 2016-01-25 15:44:15 +01:00
parent ef81bcb0b4
commit 4696e3d5a9
8 changed files with 1090 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

552
data/pseudo/bfd Normal file
View File

@ -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

2
ocaml/.gitignore vendored
View File

@ -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

12
ocaml/Gto.ml
View File

@ -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"
;;

11
ocaml/Makefile
View File

@ -14,16 +14,17 @@ 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)
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:

273
ocaml/Pseudo.ml Normal file
View File

@ -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

242
ocaml/qp_create_ezfio_from_xyz.ml
View File

@ -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,10 +135,6 @@ 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
^ "." ^ basis ^ "\" \"" ^ basis ^"\""
in
@ -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,6 +536,7 @@ 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
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
@ -339,15 +547,7 @@ let run ?o b c d m p xyz_file =
~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;
in
match Input.Ao_basis.read () with
| None -> failwith "Error in basis"
| Some x -> Input.Ao_basis.write x

24
ocaml/qp_edit.ml
View File

@ -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 ;
]

6
ocaml/qptypes_generator.ml
View File

@ -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
"
;;