From 60e374eb0365f0ba768cc9a7f09ebe2009791cb3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 19 Jan 2018 17:42:12 +0100 Subject: [PATCH] ERI tested with QP -> OK --- Basis/Basis.ml | 36 +++++---- Basis/Contracted_shell.ml | 16 +++- Basis/ERI.ml | 162 ++++++++++++++++++++++++++++---------- Makefile | 3 +- Utils/Constants.ml | 8 ++ Utils/Coordinate.ml | 52 ++++++------ Utils/Util.ml | 6 +- Utils/Zkey.ml | 25 +++--- 8 files changed, 201 insertions(+), 107 deletions(-) create mode 100644 Utils/Constants.ml diff --git a/Basis/Basis.ml b/Basis/Basis.ml index bef0162..5f9ddc1 100644 --- a/Basis/Basis.ml +++ b/Basis/Basis.ml @@ -2,18 +2,26 @@ type t = Contracted_shell.t array (** Returns an array of the basis set per atom *) let of_nuclei_and_general_basis n b = - Array.map (fun (e, center) -> - List.assoc e b - |> Array.map (fun (totAngMom, shell) -> - let expo = Array.map (fun General_basis.{exponent ; coefficient} -> - exponent) shell - and coef = Array.map (fun General_basis.{exponent ; coefficient} -> - coefficient) shell - in - Contracted_shell.create ~expo ~coef ~totAngMom ~center) - ) n - |> Array.to_list - |> Array.concat + let result = + Array.map (fun (e, center) -> + List.assoc e b + |> Array.map (fun (totAngMom, shell) -> + let expo = Array.map (fun General_basis.{exponent ; coefficient} -> + exponent) shell + and coef = Array.map (fun General_basis.{exponent ; coefficient} -> + coefficient) shell + in + Contracted_shell.create ~expo ~coef ~totAngMom ~center ~indice:0) + ) n + |> Array.to_list + |> 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 = @@ -25,8 +33,8 @@ let to_string b = ---------------- ----------------------------------------------------------------------- -Angular Coordinates (Bohr) Exponents Coefficients -Momentum X Y Z + # Angular Coordinates (Bohr) Exponents Coefficients + Momentum X Y Z ----------------------------------------------------------------------- " ^ diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml index 02d8558..389aad0 100644 --- a/Basis/Contracted_shell.ml +++ b/Basis/Contracted_shell.ml @@ -7,6 +7,8 @@ type t = { totAngMom : Angular_momentum.t; size : int; norm_coef : (int array -> float) array; + indice : int; + powers : Zkey.t array; } let size a = a.size @@ -15,6 +17,8 @@ let coef a i = a.coef.(i) let center a = a.center let totAngMom a = a.totAngMom let norm_coef a i = a.norm_coef.(i) +let indice a = a.indice +let powers a = a.powers let to_string s = @@ -22,7 +26,11 @@ let to_string s = Coordinate.to_Bohr s.center 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) ) ^ (Array.map2 (fun e c -> sprintf "%16.8e %16.8e" e c) s.expo s.coef |> Array.to_list |> String.concat (sprintf "\n%36s" " ") ) @@ -50,11 +58,13 @@ let compute_norm_coef s = ) 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 > 0); 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 { tmp with norm_coef = compute_norm_coef tmp } + diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 467422c..9d6d6d7 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -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 } -(** Write all integrals to a file *) + +(** Write all integrals to a file with the convention *) let to_file ~filename basis = let oc = open_out filename 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 - print_int !i_shift ; print_newline (); - j_shift := 1; + print_int basis.(i).Contracted_shell.indice ; print_newline (); for j=0 to i do - k_shift := 1; for k=0 to i do - l_shift := 1; for l=0 to k do + (* Compute all the integrals of the class *) let cls = contracted_class basis.(i) basis.(j) basis.(k) basis.(l) in - for i_c = 0 to zkey.(i).n - 1 do - let x = Zkey.(to_int_array Kind_3 zkey.(i).cls.(i_c)) in + for i_c = 0 to (Array.length basis.(i).Contracted_shell.powers) - 1 do + let x = Zkey.(to_int_array Kind_3 basis.(i).Contracted_shell.powers.(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 + for j_c = 0 to (Array.length basis.(j).Contracted_shell.powers) - 1 do + let x = Zkey.(to_int_array Kind_3 basis.(j).Contracted_shell.powers.(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 + for k_c = 0 to (Array.length basis.(k).Contracted_shell.powers) - 1 do + let x = Zkey.(to_int_array Kind_3 basis.(k).Contracted_shell.powers.(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 + for l_c = 0 to (Array.length basis.(l).Contracted_shell.powers) - 1 do + let x = Zkey.(to_int_array Kind_3 basis.(l).Contracted_shell.powers.(l_c)) in key_array.( 9) <- x.(0); key_array.(10) <- x.(1); key_array.(11) <- x.(2); @@ -86,25 +74,45 @@ let to_file ~filename basis = in if (abs_float value > cutoff) then 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 done done done done; - l_shift := !l_shift + zkey.(l).n done; - k_shift := !k_shift + zkey.(k).n done; - j_shift := !j_shift + zkey.(j).n done; - i_shift := !i_shift + zkey.(i).n done ; 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 oc = open_out filename in let zkey = Array.map (fun b -> @@ -159,9 +167,7 @@ let to_file ~filename basis = in if (abs_float value > cutoff) then let key = - Zkey.of_int_array Zkey.Kind_4 [| - (!i_shift+i_c);(!j_shift+j_c);(!k_shift+k_c);(!l_shift+l_c) - |] + index (!i_shift+i_c) (!j_shift+j_c) (!k_shift+k_c) (!l_shift+l_c) in result := (key, value) :: !result done @@ -177,27 +183,97 @@ let to_file ~filename basis = i_shift := !i_shift + zkey.(i).n done ; + print_endline "Computation Done"; 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 j=0 to !i_shift - 2 do - for k=0 to !i_shift - 2 do - for l=0 to !i_shift - 2 do + for i=1 to !i_shift - 1 do + for k=1 to !i_shift - 1 do + for j=1 to !i_shift - 1 do + for l=1 to !i_shift - 1 do let key = - Zkey.of_int_array Zkey.Kind_4 [| i;j;k;l |] + index i j k l in try let value = - Zmap.find result key + IntegralMap.find result key in - Printf.fprintf oc "%4d %4d %4d %4d %20.12e\n" - i j k l value + Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i j k l value with Not_found -> () done done done done; 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 + ; + *) + + diff --git a/Makefile b/Makefile index dee7889..20bc775 100644 --- a/Makefile +++ b/Makefile @@ -4,7 +4,6 @@ INCLUDE_DIRS=Nuclei,Utils,Basis LIBS= PKGS= OCAMLCFLAGS="-g" -OCAMLCFLAGS="-unsafe -noassert -safe-string" OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) -Is $(INCLUDE_DIRS) MLLFILES=$(wildcard */*.mll) $(wildcard *.mll) MLYFILES=$(wildcard */*.mly) $(wildcard *.mly) @@ -33,7 +32,7 @@ doc: qpackage.odocl $(OCAMLBUILD) $*.byte -use-ocamlfind $(PKGS) ln -s $*.byte $* -%.native: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) +%.native: $(MLFILES) $(MLIFILES) $(MLLFILES) $(MLYFILES) %.byte rm -f -- $* $(OCAMLBUILD) $*.native -use-ocamlfind $(PKGS) ln -s $*.native $* diff --git a/Utils/Constants.ml b/Utils/Constants.ml new file mode 100644 index 0000000..c2d1085 --- /dev/null +++ b/Utils/Constants.ml @@ -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 diff --git a/Utils/Coordinate.ml b/Utils/Coordinate.ml index 9c0dbe6..51cd7bf 100644 --- a/Utils/Coordinate.ml +++ b/Utils/Coordinate.ml @@ -2,9 +2,7 @@ type t = | Bohr of float array | Angstrom of float array -(** Bohr radius, taken from https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0 *) -let a0 = - 0.529_177_210_67 +let a0 = Constants.a0 let zero = Bohr [| 0. ; 0. ; 0. |] @@ -24,21 +22,10 @@ let to_string y = | Angstrom x -> (result x) ^ " Angstrom" -let to_float_array = function +let extract_float_array = function | Bohr 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 *) let (|.) s a = @@ -49,10 +36,10 @@ let (|.) s a = let to_Angstrom = function | 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 - | Angstrom a -> Bohr (1./.a0 |. Angstrom a |> to_float_array) + | Angstrom a -> Bohr (1./.a0 |. Angstrom a |> extract_float_array) | Bohr a -> Bohr a let (|-), (|+) = @@ -76,22 +63,29 @@ let (|-), (|+) = ) -let dot = - let rec op f p q = - match (p, q) with - | (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) +let dot p q = + let f = function + | Bohr [|x;y;z|], Bohr [|x';y';z'|] -> x*.x' +. y*.y' +. z*.z' + | _ -> assert false in - op (fun a b -> - match a,b with - | [|x;y;z|], [|x';y';z'|] -> x*.x' +. y*.y' +. z*.z' - | _ -> assert false - ) + f (to_Bohr p, to_Bohr q) let norm 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") + diff --git a/Utils/Util.ml b/Utils/Util.ml index 5a666ff..8ebca5b 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -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 diff --git a/Utils/Zkey.ml b/Utils/Zkey.ml index a3820b6..28e0d25 100644 --- a/Utils/Zkey.ml +++ b/Utils/Zkey.ml @@ -3,14 +3,15 @@ include Z type kind = -| Kind_2 | Kind_3 -| Kind_4 | Kind_6 | 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 (<|) x 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) in 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_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) |> Z.of_int64 | 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) |> Z.of_int64 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 *) let to_int_array ~kind a = 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 ; Z.to_int @@ Z.extract a 10 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 ; Z.to_int @@ Z.extract a 40 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 10 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 = "< " ^ ( Z.to_string a ) ^ " | " ^ (