LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-11-06 22:23:42 +01:00
Code Releases Activity

ERI tested with QP -> OK

Browse Source
This commit is contained in:
Anthony Scemama 2018-01-19 17:42:12 +01:00
parent 73c8e48731
commit 60e374eb03
8 changed files with 201 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Basis/Basis.ml
View File

@ -2,6 +2,7 @@ type t = Contracted_shell.t array
(** Returns an array of the basis set per atom *) (** Returns an array of the basis set per atom *)
let of_nuclei_and_general_basis n b = let of_nuclei_and_general_basis n b =
let result =
Array.map (fun (e, center) -> Array.map (fun (e, center) ->
List.assoc e b List.assoc e b
|> Array.map (fun (totAngMom, shell) -> |> Array.map (fun (totAngMom, shell) ->
@ -10,10 +11,17 @@ let of_nuclei_and_general_basis n b =
and coef = Array.map (fun General_basis.{exponent ; coefficient} -> and coef = Array.map (fun General_basis.{exponent ; coefficient} ->
coefficient) shell coefficient) shell
in in
Contracted_shell.create ~expo ~coef ~totAngMom ~center) Contracted_shell.create ~expo ~coef ~totAngMom ~center ~indice:0)
) n ) n
|> Array.to_list |> Array.to_list
|> Array.concat |> Array.concat
in
Array.iteri (fun i x ->
if (i > 0) then
result.(i) <- {x with Contracted_shell.indice= (
result.(i-1).Contracted_shell.indice + (Array.length result.(i-1).Contracted_shell.powers)) }
) result ;
result
let to_string b = let to_string b =
@ -25,7 +33,7 @@ let to_string b =
---------------- ----------------
----------------------------------------------------------------------- -----------------------------------------------------------------------
Angular Coordinates (Bohr) Exponents Coefficients # Angular Coordinates (Bohr) Exponents Coefficients
Momentum X Y Z Momentum X Y Z
----------------------------------------------------------------------- -----------------------------------------------------------------------
" "

16
Basis/Contracted_shell.ml
View File

@ -7,6 +7,8 @@ type t = {
totAngMom : Angular_momentum.t; totAngMom : Angular_momentum.t;
size : int; size : int;
norm_coef : (int array -> float) array; norm_coef : (int array -> float) array;
indice : int;
powers : Zkey.t array;
} }
let size a = a.size let size a = a.size
@ -15,6 +17,8 @@ let coef a i = a.coef.(i)
let center a = a.center let center a = a.center
let totAngMom a = a.totAngMom let totAngMom a = a.totAngMom
let norm_coef a i = a.norm_coef.(i) let norm_coef a i = a.norm_coef.(i)
let indice a = a.indice
let powers a = a.powers
let to_string s = let to_string s =
@ -22,7 +26,11 @@ let to_string s =
Coordinate.to_Bohr s.center Coordinate.to_Bohr s.center
in in
let open Printf in let open Printf in
( sprintf "%2s %10.6f %10.6f %10.6f " (Angular_momentum.to_string s.totAngMom) (match s.totAngMom with
| Angular_momentum.S -> sprintf "%3d " (s.indice+1)
| _ -> sprintf "%3d-%-3d" (s.indice+1) (s.indice+(Array.length s.powers))
) ^
( sprintf "%1s %8.3f %8.3f %8.3f " (Angular_momentum.to_string s.totAngMom)
(Coordinate.x coord) (Coordinate.y coord) (Coordinate.z coord) ) ^ (Coordinate.x coord) (Coordinate.y coord) (Coordinate.z coord) ) ^
(Array.map2 (fun e c -> sprintf "%16.8e %16.8e" e c) s.expo s.coef (Array.map2 (fun e c -> sprintf "%16.8e %16.8e" e c) s.expo s.coef
|> Array.to_list |> String.concat (sprintf "\n%36s" " ") ) |> Array.to_list |> String.concat (sprintf "\n%36s" " ") )
@ -50,11 +58,13 @@ let compute_norm_coef s =
) s.expo ) s.expo
let create ~expo ~coef ~center ~totAngMom = let create ~indice ~expo ~coef ~center ~totAngMom =
assert (Array.length expo = Array.length coef); assert (Array.length expo = Array.length coef);
assert (Array.length expo > 0); assert (Array.length expo > 0);
let tmp = let tmp =
{ expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef = [||]} { indice ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef = [||];
powers = Angular_momentum.zkey_array (Kind_1 totAngMom) }
in in
{ tmp with norm_coef = compute_norm_coef tmp } { tmp with norm_coef = compute_norm_coef tmp }

160
Basis/ERI.ml
View File

@ -31,50 +31,38 @@ let contracted_class shell_a shell_b shell_c shell_d : float Zmap.t =
type n_cls = { n : int ; cls : Z.t array } type n_cls = { n : int ; cls : Z.t array }
(** Write all integrals to a file *)
(** Write all integrals to a file with the <ij|kl> convention *)
let to_file ~filename basis = let to_file ~filename basis =
let oc = open_out filename in let oc = open_out filename in
let key_array = Array.make 12 0 in let key_array = Array.make 12 0 in
let zkey = Array.map (fun b ->
let result =
Angular_momentum.(zkey_array (Kind_1 b.Contracted_shell.totAngMom))
in
{ n=Array.length result ; cls=result }
) basis
in
let i_shift = ref 1 in
let j_shift = ref 1 in
let k_shift = ref 1 in
let l_shift = ref 1 in
for i=0 to (Array.length basis) - 1 do for i=0 to (Array.length basis) - 1 do
print_int !i_shift ; print_newline (); print_int basis.(i).Contracted_shell.indice ; print_newline ();
j_shift := 1;
for j=0 to i do for j=0 to i do
k_shift := 1;
for k=0 to i do for k=0 to i do
l_shift := 1;
for l=0 to k do for l=0 to k do
(* Compute all the integrals of the class *)
let cls = let cls =
contracted_class basis.(i) basis.(j) basis.(k) basis.(l) contracted_class basis.(i) basis.(j) basis.(k) basis.(l)
in in
for i_c = 0 to zkey.(i).n - 1 do for i_c = 0 to (Array.length basis.(i).Contracted_shell.powers) - 1 do
let x = Zkey.(to_int_array Kind_3 zkey.(i).cls.(i_c)) in let x = Zkey.(to_int_array Kind_3 basis.(i).Contracted_shell.powers.(i_c)) in
key_array.(0) <- x.(0); key_array.(0) <- x.(0);
key_array.(1) <- x.(1); key_array.(1) <- x.(1);
key_array.(2) <- x.(2); key_array.(2) <- x.(2);
for j_c = 0 to zkey.(j).n - 1 do for j_c = 0 to (Array.length basis.(j).Contracted_shell.powers) - 1 do
let x = Zkey.(to_int_array Kind_3 zkey.(j).cls.(j_c)) in let x = Zkey.(to_int_array Kind_3 basis.(j).Contracted_shell.powers.(j_c)) in
key_array.(3) <- x.(0); key_array.(3) <- x.(0);
key_array.(4) <- x.(1); key_array.(4) <- x.(1);
key_array.(5) <- x.(2); key_array.(5) <- x.(2);
for k_c = 0 to zkey.(k).n - 1 do for k_c = 0 to (Array.length basis.(k).Contracted_shell.powers) - 1 do
let x = Zkey.(to_int_array Kind_3 zkey.(k).cls.(k_c)) in let x = Zkey.(to_int_array Kind_3 basis.(k).Contracted_shell.powers.(k_c)) in
key_array.(6) <- x.(0); key_array.(6) <- x.(0);
key_array.(7) <- x.(1); key_array.(7) <- x.(1);
key_array.(8) <- x.(2); key_array.(8) <- x.(2);
for l_c = 0 to zkey.(l).n - 1 do for l_c = 0 to (Array.length basis.(l).Contracted_shell.powers) - 1 do
let x = Zkey.(to_int_array Kind_3 zkey.(l).cls.(l_c)) in let x = Zkey.(to_int_array Kind_3 basis.(l).Contracted_shell.powers.(l_c)) in
key_array.( 9) <- x.(0); key_array.( 9) <- x.(0);
key_array.(10) <- x.(1); key_array.(10) <- x.(1);
key_array.(11) <- x.(2); key_array.(11) <- x.(2);
@ -86,25 +74,45 @@ let to_file ~filename basis =
in in
if (abs_float value > cutoff) then if (abs_float value > cutoff) then
Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n"
(!i_shift+i_c) (!j_shift+j_c) (!k_shift+k_c) (!l_shift+l_c) (basis.(i).Contracted_shell.indice + i_c + 1)
(basis.(k).Contracted_shell.indice + k_c + 1)
(basis.(j).Contracted_shell.indice + j_c + 1)
(basis.(l).Contracted_shell.indice + l_c + 1)
value value
done done
done done
done done
done; done;
l_shift := !l_shift + zkey.(l).n
done; done;
k_shift := !k_shift + zkey.(k).n
done; done;
j_shift := !j_shift + zkey.(j).n
done; done;
i_shift := !i_shift + zkey.(i).n
done done
; ;
close_out oc close_out oc
(* (*
module IntegralMap = Zmap
let index i j k l =
Zkey.of_int_array Zkey.Kind_4 [| i;j;k;l |]
module Key = struct
type t=int
let equal (x:int) (y:int) = x = y
let hash (x:int) = x
end
module IntegralMap = Hashtbl.Make(Key)
let index i j k l =
let f i k =
let (p,r) =
if i <= k then (i,k) else (k,i)
in p+ (r*r-r)/2
in
let p = f i k and q = f j l in
f p q
let to_file ~filename basis = let to_file ~filename basis =
let oc = open_out filename in let oc = open_out filename in
let zkey = Array.map (fun b -> let zkey = Array.map (fun b ->
@ -159,9 +167,7 @@ let to_file ~filename basis =
in in
if (abs_float value > cutoff) then if (abs_float value > cutoff) then
let key = let key =
Zkey.of_int_array Zkey.Kind_4 [| index (!i_shift+i_c) (!j_shift+j_c) (!k_shift+k_c) (!l_shift+l_c)
(!i_shift+i_c);(!j_shift+j_c);(!k_shift+k_c);(!l_shift+l_c)
|]
in in
result := (key, value) :: !result result := (key, value) :: !result
done done
@ -177,27 +183,97 @@ let to_file ~filename basis =
i_shift := !i_shift + zkey.(i).n i_shift := !i_shift + zkey.(i).n
done done
; ;
print_endline "Computation Done";
let result = Array.of_list !result in let result = Array.of_list !result in
let result = Zmap.create (Array.length result) in let result =
let a = IntegralMap.create (Array.length result) in
Array.iter (fun (k,v) -> IntegralMap.add a k v) result;
a
in
print_endline "Map formed";
for i=0 to !i_shift - 2 do for i=1 to !i_shift - 1 do
for j=0 to !i_shift - 2 do for k=1 to !i_shift - 1 do
for k=0 to !i_shift - 2 do for j=1 to !i_shift - 1 do
for l=0 to !i_shift - 2 do for l=1 to !i_shift - 1 do
let key = let key =
Zkey.of_int_array Zkey.Kind_4 [| i;j;k;l |] index i j k l
in in
try try
let value = let value =
Zmap.find result key IntegralMap.find result key
in in
Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i j k l value
i j k l value
with Not_found -> () with Not_found -> ()
done done
done done
done done
done; done;
close_out oc close_out oc
let xto_file ~filename basis =
let zkey = Array.map (fun b ->
let result =
Angular_momentum.(zkey_array (Kind_1 b.Contracted_shell.totAngMom))
in
{ n=Array.length result ; cls=result }
) basis
in
let key_array = Array.make 12 0 in
let result = ref [] in
let (i,j,k,l) = (1,1,1,18) in
let (i,j,k,l) = (i-1,j-1,k-1,l-1) in
basis.(i) |> Contracted_shell.to_string |> print_endline;
basis.(j) |> Contracted_shell.to_string |> print_endline;
basis.(k) |> Contracted_shell.to_string |> print_endline;
basis.(l) |> Contracted_shell.to_string |> print_endline;
let bi, bj, bk, bl =
basis.(i), basis.(j), basis.(k), basis.(l)
in
let cls =
(*contracted_class basis.(i) basis.(j) basis.(k) basis.(l) *)
contracted_class bi bj bk bl
in
Zmap.iter (fun k v -> Printf.printf "%50s %20e\n" Zkey.(to_string Kind_12 k) v) cls;
for i_c = 0 to zkey.(i).n - 1 do
let x = Zkey.(to_int_array Kind_3 zkey.(i).cls.(i_c)) in
key_array.(0) <- x.(0);
key_array.(1) <- x.(1);
key_array.(2) <- x.(2);
for j_c = 0 to zkey.(j).n - 1 do
let x = Zkey.(to_int_array Kind_3 zkey.(j).cls.(j_c)) in
key_array.(3) <- x.(0);
key_array.(4) <- x.(1);
key_array.(5) <- x.(2);
for k_c = 0 to zkey.(k).n - 1 do
let x = Zkey.(to_int_array Kind_3 zkey.(k).cls.(k_c)) in
key_array.(6) <- x.(0);
key_array.(7) <- x.(1);
key_array.(8) <- x.(2);
for l_c = 0 to zkey.(l).n - 1 do
let x = Zkey.(to_int_array Kind_3 zkey.(l).cls.(l_c)) in
key_array.( 9) <- x.(0);
key_array.(10) <- x.(1);
key_array.(11) <- x.(2);
let key =
Zkey.(of_int_array Kind_12 key_array)
in
let value =
Zmap.find cls key
in
if (abs_float value > cutoff) then
result := (key, value) :: !result
done
done
done
done;
List.iter (fun (k,v) -> Printf.printf "%60s %e\n" Zkey.(to_string Kind_12 k) v) !result
;
*) *)

3
Makefile
View File

@ -4,7 +4,6 @@ INCLUDE_DIRS=Nuclei,Utils,Basis
LIBS= LIBS=
PKGS= PKGS=
OCAMLCFLAGS="-g" OCAMLCFLAGS="-g"
OCAMLCFLAGS="-unsafe -noassert -safe-string"
OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) -Is $(INCLUDE_DIRS) OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) -Is $(INCLUDE_DIRS)
MLLFILES=$(wildcard */*.mll) $(wildcard *.mll) MLLFILES=$(wildcard */*.mll) $(wildcard *.mll)
MLYFILES=$(wildcard */*.mly) $(wildcard *.mly) MLYFILES=$(wildcard */*.mly) $(wildcard *.mly)
@ -33,7 +32,7 @@ doc: qpackage.odocl
$(OCAMLBUILD) $*.byte -use-ocamlfind $(PKGS) $(OCAMLBUILD) $*.byte -use-ocamlfind $(PKGS)
ln -s $*.byte $* ln -s $*.byte $*
%.native: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) %.native: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) %.byte
rm -f -- $* rm -f -- $*
$(OCAMLBUILD) $*.native -use-ocamlfind $(PKGS) $(OCAMLBUILD) $*.native -use-ocamlfind $(PKGS)
ln -s $*.native $* ln -s $*.native $*

8
Utils/Constants.ml Normal file
View File

@ -0,0 +1,8 @@
let cutoff = 1.e-16
(** Constants *)
let pi = acos (-1.)
let pi_inv = 1. /. pi
let two_over_sq_pi = 2. /. (sqrt pi)
let a0 = 0.529_177_210_67

52
Utils/Coordinate.ml
View File

@ -2,9 +2,7 @@ type t =
| Bohr of float array | Bohr of float array
| Angstrom of float array | Angstrom of float array
(** Bohr radius, taken from https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0 *) let a0 = Constants.a0
let a0 =
0.529_177_210_67
let zero = Bohr [| 0. ; 0. ; 0. |] let zero = Bohr [| 0. ; 0. ; 0. |]
@ -24,21 +22,10 @@ let to_string y =
| Angstrom x -> (result x) ^ " Angstrom" | Angstrom x -> (result x) ^ " Angstrom"
let to_float_array = function let extract_float_array = function
| Bohr a | Bohr a
| Angstrom a -> a | Angstrom a -> a
let x a = (to_float_array a).(0)
let y a = (to_float_array a).(1)
let z a = (to_float_array a).(2)
let coord a = function
| 0 -> (to_float_array a).(0)
| 1 -> (to_float_array a).(1)
| 2 -> (to_float_array a).(2)
| _ -> raise (Invalid_argument "Coordinate")
(** Linear algebra *) (** Linear algebra *)
let (|.) s a = let (|.) s a =
@ -49,10 +36,10 @@ let (|.) s a =
let to_Angstrom = function let to_Angstrom = function
| Angstrom a -> Angstrom a | Angstrom a -> Angstrom a
| Bohr a -> Angstrom (a0 |. Bohr a |> to_float_array) | Bohr a -> Angstrom (a0 |. Bohr a |> extract_float_array)
let to_Bohr = function let to_Bohr = function
| Angstrom a -> Bohr (1./.a0 |. Angstrom a |> to_float_array) | Angstrom a -> Bohr (1./.a0 |. Angstrom a |> extract_float_array)
| Bohr a -> Bohr a | Bohr a -> Bohr a
let (|-), (|+) = let (|-), (|+) =
@ -76,22 +63,29 @@ let (|-), (|+) =
) )
let dot = let dot p q =
let rec op f p q = let f = function
match (p, q) with | Bohr [|x;y;z|], Bohr [|x';y';z'|] -> x*.x' +. y*.y' +. z*.z'
| (Angstrom a, Angstrom b)
| (Bohr a, Bohr b) -> f a b
| (Angstrom a, Bohr b) -> op f (to_Bohr p) q
| (Bohr a, Angstrom b) -> op f p (to_Bohr q)
in
op (fun a b ->
match a,b with
| [|x;y;z|], [|x';y';z'|] -> x*.x' +. y*.y' +. z*.z'
| _ -> assert false | _ -> assert false
) in
f (to_Bohr p, to_Bohr q)
let norm u = let norm u =
sqrt @@ dot u u sqrt @@ dot u u
let rec to_float_array a =
to_Bohr a |> extract_float_array
let x a = (extract_float_array @@ to_Bohr a).(0)
let y a = (extract_float_array @@ to_Bohr a).(1)
let z a = (extract_float_array @@ to_Bohr a).(2)
let coord a = function
| 0 -> (extract_float_array @@ to_Bohr a).(0)
| 1 -> (extract_float_array @@ to_Bohr a).(1)
| 2 -> (extract_float_array @@ to_Bohr a).(2)
| _ -> raise (Invalid_argument "Coordinate")

6
Utils/Util.ml
View File

@ -1,9 +1,5 @@
let cutoff = 1.e-50 include Constants
(** Constants *)
let pi = acos (-1.)
let pi_inv = 1. /. pi
let two_over_sq_pi = 2. /. (sqrt pi)
let factmax = 150 let factmax = 150

25
Utils/Zkey.ml
View File

@ -3,14 +3,15 @@
include Z include Z
type kind = type kind =
| Kind_2
| Kind_3 | Kind_3
| Kind_4
| Kind_6 | Kind_6
| Kind_12 | Kind_12
| Kind_4
| Kind_2
| Kind_1
(** Build a Zkey from an array or 2, 3, 4, 6, or 12 integers *) (** Build a Zkey from an array or 1, 2, 3, 4, 6, or 12 integers *)
let of_int_array ~kind a = let of_int_array ~kind a =
let (<|) x a = let (<|) x a =
Z.logor (Z.shift_left x 64) a Z.logor (Z.shift_left x 64) a
@ -22,9 +23,7 @@ let of_int_array ~kind a =
Int64.logor (Int64.shift_left x 16) (Int64.of_int a) Int64.logor (Int64.shift_left x 16) (Int64.of_int a)
in in
match kind with match kind with
| Kind_2 -> (Int64.of_int a.(0)) <+ a.(1) |> Z.of_int64
| Kind_3 -> (Int64.of_int a.(0)) << a.(1) << a.(2) |> Z.of_int64 | Kind_3 -> (Int64.of_int a.(0)) << a.(1) << a.(2) |> Z.of_int64
| Kind_4 -> (Int64.of_int a.(0)) <+ a.(1) <+ a.(2) <+ a.(3) |> Z.of_int64
| Kind_6 -> (Int64.of_int a.(0)) << a.(1) << a.(2) << a.(3) << a.(4) << a.(5) | Kind_6 -> (Int64.of_int a.(0)) << a.(1) << a.(2) << a.(3) << a.(4) << a.(5)
|> Z.of_int64 |> Z.of_int64
| Kind_12 -> | Kind_12 ->
@ -35,20 +34,17 @@ let of_int_array ~kind a =
(Int64.of_int a.(6)) << a.(7) << a.(8) << a.(9) << a.(10) << a.(11) (Int64.of_int a.(6)) << a.(7) << a.(8) << a.(9) << a.(10) << a.(11)
|> Z.of_int64 |> Z.of_int64
in a <| b in a <| b
| Kind_4 -> (Int64.of_int a.(0)) <+ a.(1) <+ a.(2) <+ a.(3) |> Z.of_int64
| Kind_2 -> (Int64.of_int a.(0)) <+ a.(1) |> Z.of_int64
| Kind_1 -> Z.of_int a.(0)
(** Transform the Zkey into an int array *) (** Transform the Zkey into an int array *)
let to_int_array ~kind a = let to_int_array ~kind a =
match kind with match kind with
| Kind_2 -> [| Z.to_int @@ Z.extract a 16 16 ;
Z.to_int @@ Z.extract a 0 16 |]
| Kind_3 -> [| Z.to_int @@ Z.extract a 20 10 ; | Kind_3 -> [| Z.to_int @@ Z.extract a 20 10 ;
Z.to_int @@ Z.extract a 10 10 ; Z.to_int @@ Z.extract a 10 10 ;
Z.to_int @@ Z.extract a 0 10 |] Z.to_int @@ Z.extract a 0 10 |]
| Kind_4 -> [| Z.to_int @@ Z.extract a 48 16 ;
Z.to_int @@ Z.extract a 32 16 ;
Z.to_int @@ Z.extract a 16 16 ;
Z.to_int @@ Z.extract a 0 16 |]
| Kind_6 -> [| Z.to_int @@ Z.extract a 50 10 ; | Kind_6 -> [| Z.to_int @@ Z.extract a 50 10 ;
Z.to_int @@ Z.extract a 40 10 ; Z.to_int @@ Z.extract a 40 10 ;
Z.to_int @@ Z.extract a 30 10 ; Z.to_int @@ Z.extract a 30 10 ;
@ -67,6 +63,13 @@ let to_int_array ~kind a =
Z.to_int @@ Z.extract a 20 10 ; Z.to_int @@ Z.extract a 20 10 ;
Z.to_int @@ Z.extract a 10 10 ; Z.to_int @@ Z.extract a 10 10 ;
Z.to_int @@ Z.extract a 0 10 |] Z.to_int @@ Z.extract a 0 10 |]
| Kind_4 -> [| Z.to_int @@ Z.extract a 48 16 ;
Z.to_int @@ Z.extract a 32 16 ;
Z.to_int @@ Z.extract a 16 16 ;
Z.to_int @@ Z.extract a 0 16 |]
| Kind_2 -> [| Z.to_int @@ Z.extract a 16 16 ;
Z.to_int @@ Z.extract a 0 16 |]
| Kind_1 -> [| Z.to_int a |]
let to_string ~kind a = let to_string ~kind a =
"< " ^ ( Z.to_string a ) ^ " | " ^ ( "< " ^ ( Z.to_string a ) ^ " | " ^ (