From 6437c5503acc890b0f76d1b506f8ef44de077884 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 17 Jan 2018 18:19:38 +0100 Subject: [PATCH] Working on AO shells --- Basis/Contracted_shell.ml | 44 ++++++++++++++++ Basis/Shell_pair.ml | 103 ++++++++++++++++++++++++++++++++++++++ Utils/Angular_momentum.ml | 65 ++++++++++++++++++++++++ Utils/Coordinate.ml | 55 ++++++++++++++++++++ Utils/Coordinate.mli | 16 ++++++ Utils/Zkey.ml | 96 +++++++++++++++++++++++++++++++++++ Utils/Zmap.ml | 4 ++ 7 files changed, 383 insertions(+) create mode 100644 Basis/Contracted_shell.ml create mode 100644 Basis/Shell_pair.ml create mode 100644 Utils/Coordinate.ml create mode 100644 Utils/Coordinate.mli create mode 100644 Utils/Zkey.ml create mode 100644 Utils/Zmap.ml diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml new file mode 100644 index 0000000..3138d5c --- /dev/null +++ b/Basis/Contracted_shell.ml @@ -0,0 +1,44 @@ +open Util + +type t = { + expo : float array; + coef : float array; + center : Coordinate.t; + totAngMom : Angular_momentum.t; + size : int; + norm_coef : (int array -> float) array; +} + + + +(** Normalization coefficient of contracted function i, which depends on the + exponent and the angular momentum. Two conventions can be chosen : a single + normalisation factor for all functions of the class, or a coefficient which + depends on the powers of x,y and z. +*) +let compute_norm_coef s = + let atot = + Angular_momentum.to_int s.totAngMom + in + Array.mapi (fun i alpha -> + let c = + ((alpha +. alpha) *. pi_inv)**(1.5) *. (pow (4. *. alpha) atot) + in + let result a = + let dfa = Array.map (fun j -> + fact (j+j) /. ( float_of_int (1 lsl j) *. fact j) + ) a + in sqrt (c /. (dfa.(0) *.dfa.(1) *. dfa.(2))) + in + result + ) s.expo + + +let create ~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 = [||]} + in + { tmp with norm_coef = compute_norm_coef tmp } + diff --git a/Basis/Shell_pair.ml b/Basis/Shell_pair.ml new file mode 100644 index 0000000..d2fb9e5 --- /dev/null +++ b/Basis/Shell_pair.ml @@ -0,0 +1,103 @@ +open Util + +type t = { + expo : float; + expo_inv : float; + center_ab: float array; + center : float array; + norm_sq : float; + norm : float; + coef : float; + norm_fun : int array -> int array -> float; + i : int; + j : int; +} + +exception Null_contribution + +let create_array ?(cutoff=0.) p_a p_b = + let log_cutoff = + if (cutoff = 0.) then infinity + else -. (log cutoff) + in + + let x_a = Coordinate.x p_a.Contracted_shell.center + and y_a = Coordinate.y p_a.Contracted_shell.center + and z_a = Coordinate.z p_a.Contracted_shell.center + and x_b = Coordinate.x p_b.Contracted_shell.center + and y_b = Coordinate.y p_b.Contracted_shell.center + and z_b = Coordinate.z p_b.Contracted_shell.center + in + (* + match p_a.Contracted_shell.center, p_b.Contracted_shell.center with + | [|x_a; y_a; z_a|], [|x_b; y_b; z_b|] -> + *) + let center_ab = + Coordinate.(p_a.Contracted_shell.center |- p_b.Contracted_shell.center) + in + let norm_sq = + Coordinate.dot center_ab center_ab + in + Array.init p_a.Contracted_shell.size (fun i -> + let p_a_expo_center = + [| p_a.Contracted_shell.expo.(i) *. x_a ; p_a.Contracted_shell.expo.(i) *. y_a ; p_a.Contracted_shell.expo.(i) *. z_a |] + in + + Array.init p_b.Contracted_shell.size (fun j -> + try + let f1 = + p_a.Contracted_shell.norm_coef.(i) + in + let f2 = + p_b.Contracted_shell.norm_coef.(j) + in + let norm_fun a b = + f1 a *. f2 b + in + let norm = + norm_fun + [| Angular_momentum.to_int p_a.Contracted_shell.totAngMom ; 0 ; 0 |] + [| Angular_momentum.to_int p_b.Contracted_shell.totAngMom ; 0 ; 0 |] + in + if (norm < cutoff) then + raise Null_contribution; + let p_b_expo_center = + [| p_b.Contracted_shell.expo.(j) *. x_b ; p_b.Contracted_shell.expo.(j) *. y_b ; p_b.Contracted_shell.expo.(j) *. z_b |] + in + let expo = p_a.Contracted_shell.expo.(i) +. p_b.Contracted_shell.expo.(j) in + let expo_inv = 1. /. expo in + let center = + [| (p_a_expo_center.(0) +. p_b_expo_center.(0)) *. expo_inv; + (p_a_expo_center.(1) +. p_b_expo_center.(1)) *. expo_inv; + (p_a_expo_center.(2) +. p_b_expo_center.(2)) *. expo_inv |] + in + let argexpo = + p_a.Contracted_shell.expo.(i) *. p_b.Contracted_shell.expo.(j) + *. norm_sq *. expo_inv + in + if (argexpo > log_cutoff) then + raise Null_contribution; + let g = + (pi *. expo_inv)**(1.5) *. exp(-. argexpo) + in + let norm_inv = 1./.norm in + let norm_fun a b = + norm_inv *. norm_fun a b + in + let coef = + norm *. p_a.Contracted_shell.coef.(i) *. p_b.Contracted_shell.coef.(j) *. g + in + if (abs_float coef < cutoff) then + raise Null_contribution; + Some { i ; j ; norm_fun ; norm ; coef ; expo ; expo_inv ; center ; center_ab=(Coordinate.to_float_array center_ab) ; norm_sq } + with + | Null_contribution -> None + ) + ) + |> Array.to_list + |> Array.concat + |> Array.to_list + |> List.filter (function Some _ -> true | None -> false) + |> List.map (function Some x -> x | None -> assert false) + |> Array.of_list + diff --git a/Utils/Angular_momentum.ml b/Utils/Angular_momentum.ml index 6538d98..bd48173 100644 --- a/Utils/Angular_momentum.ml +++ b/Utils/Angular_momentum.ml @@ -52,3 +52,68 @@ let of_int = function +type kind = +| Kind_2 of (t*t) +| Kind_4 of (t*t*t*t) + + +(** Returns an array of Zkeys corresponding to all possible angular momenta *) +let zkey_array a = + let keys_1d l = + let create_z (x,y,_) = + (x,y,l-(x+y)) + in + let rec create_y accu xyz = + let (x,y,z) = xyz in + match y with + | 0 -> (create_z xyz)::accu + | i -> let ynew = y-1 in + create_y ( (create_z xyz)::accu) (x,ynew,z) + in + let rec create_x accu xyz = + let (x,y,z) = xyz in + match x with + | 0 -> (create_y [] xyz)@accu + | i -> let xnew = x-1 in + let ynew = l-xnew in + create_x ((create_y [] xyz)@accu) (xnew, ynew, z) + in + create_x [] (l,0,0) + |> List.rev + in + + begin + match a with + | Kind_2 (l1, l2) -> + + let a = Array.init 6 (fun _ -> 0) in + List.map (fun (cx,cy,cz) -> + a.(0) <- cx ; a.(1) <- cy ; a.(2) <- cz ; + List.map (fun (dx,dy,dz) -> + a.(3) <- dx ; a.(4) <- dy ; a.(5) <- dz ; + Zkey.(of_int_array Kind_6 a) + ) (keys_1d @@ to_int l1) + ) (keys_1d @@ to_int l2) + + | Kind_4 (l1, l2, l3, l4) -> + + let a = Array.init 12 (fun _ -> 0) in + List.map (fun (ax,ay,az) -> + a.(0) <- ax ; a.(1) <- ay ; a.(2) <- az ; + List.map (fun (bx,by,bz) -> + a.(3) <- bx ; a.(4) <- by ; a.(5) <- bz ; + List.map (fun (cx,cy,cz) -> + a.(6) <- cx ; a.(7) <- cy ; a.(8) <- cz ; + List.map (fun (dx,dy,dz) -> + a.(9) <- dx ; a.(10) <- dy ; a.(11) <- dz ; + Zkey.(of_int_array Kind_12 a) + ) (keys_1d @@ to_int l4) + ) (keys_1d @@ to_int l3) + |> List.concat + ) (keys_1d @@ to_int l2) + |> List.concat + ) (keys_1d @@ to_int l1) + end + |> List.concat + |> Array.of_list + diff --git a/Utils/Coordinate.ml b/Utils/Coordinate.ml new file mode 100644 index 0000000..82bd610 --- /dev/null +++ b/Utils/Coordinate.ml @@ -0,0 +1,55 @@ +type t = float array + +let zero = + [| 0. ; 0. ; 0. |] + +let of_float_triplet (x,y,z) = + [|x;y;z|] + +let of_3_floats x y z = + [|x;y;z|] + +let to_string x = + (string_of_float x.(0))^" "^(string_of_float x.(1))^" "^(string_of_float x.(2)) + +let x a = a.(0) +let y a = a.(1) +let z a = a.(2) + +let coord a = function +| 0 -> a.(0) +| 1 -> a.(1) +| 2 -> a.(2) +| _ -> raise (Invalid_argument "Coordinate") + + +let to_float_array a = a + +(** Linear algebra *) +let (|-) a b = + match a,b with + | [|x;y;z|], [|x';y';z'|] -> [| x-.x'; y-.y'; z-.z' |] + | _ -> assert false + + +let (|+) a b = + match a,b with + | [|x;y;z|], [|x';y';z'|] -> [| x+.x'; y+.y'; z+.z' |] + | _ -> assert false + + +let (|.) s a = + match a with + | [|x;y;z|] -> [| s*.x; s*.y; s*.z |] + | _ -> assert false + + +let dot a b = + match a,b with + | [|x;y;z|], [|x';y';z'|] -> x*.x' +. y*.y' +. z*.z' + | _ -> assert false + + +let norm u = + sqrt @@ dot u u + diff --git a/Utils/Coordinate.mli b/Utils/Coordinate.mli new file mode 100644 index 0000000..e6bd4bb --- /dev/null +++ b/Utils/Coordinate.mli @@ -0,0 +1,16 @@ +type t +val zero : t +val of_float_triplet : float * float * float -> t +val of_3_floats : float -> float -> float -> t +val to_string : t -> string +val x : t -> float +val y : t -> float +val z : t -> float +val coord : t -> int -> float +val to_float_array : t -> float array +val (|-) : t -> t -> t +val (|+) : t -> t -> t +val (|.) : float -> t -> t +val dot : t -> t -> float +val norm : t -> float + diff --git a/Utils/Zkey.ml b/Utils/Zkey.ml new file mode 100644 index 0000000..a3820b6 --- /dev/null +++ b/Utils/Zkey.ml @@ -0,0 +1,96 @@ +(** Key for hastables that contain tuples of integers encoded in a Zarith integer *) + +include Z + +type kind = +| Kind_2 +| Kind_3 +| Kind_4 +| Kind_6 +| Kind_12 + + +(** Build a Zkey from an array or 2, 3, 4, 6, or 12 integers *) +let of_int_array ~kind a = + let (<|) x a = + Z.logor (Z.shift_left x 64) a + in + let (<<) x a = + Int64.logor (Int64.shift_left x 10) (Int64.of_int a) + in + let (<+) x 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 -> + let a = + (Int64.of_int a.(0)) << a.(1) << a.(2) << a.(3) << a.(4) << a.(5) + |> Z.of_int64 + and b = + (Int64.of_int a.(6)) << a.(7) << a.(8) << a.(9) << a.(10) << a.(11) + |> Z.of_int64 + in a <| b + + +(** 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 ; + 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_12 -> [| Z.to_int @@ Z.extract a 114 10 ; + Z.to_int @@ Z.extract a 104 10 ; + Z.to_int @@ Z.extract a 94 10 ; + Z.to_int @@ Z.extract a 84 10 ; + Z.to_int @@ Z.extract a 74 10 ; + Z.to_int @@ Z.extract a 64 10 ; + Z.to_int @@ Z.extract a 50 10 ; + Z.to_int @@ Z.extract a 40 10 ; + Z.to_int @@ Z.extract a 30 10 ; + Z.to_int @@ Z.extract a 20 10 ; + Z.to_int @@ Z.extract a 10 10 ; + Z.to_int @@ Z.extract a 0 10 |] + +let to_string ~kind a = + "< " ^ ( Z.to_string a ) ^ " | " ^ ( + to_int_array kind a + |> Array.map string_of_int + |> Array.to_list + |> String.concat ", " + ) ^ " >" + + +(* +let debug () = + let k2 = of_int_array Kind_2 [| 1 ; 2 |] + and k3 = of_int_array Kind_3 [| 1 ; 2 ; 3 |] + and k4 = of_int_array Kind_4 [| 1 ; 2 ; 3; 4 |] + and k6 = of_int_array Kind_6 [| 1 ; 2 ; 3; 4 ; 5; 6|] + and k12 = of_int_array Kind_12 [| 1 ; 2 ; 3; 4 ; 5; 6 ; 7 ; 8 ; 9 ; 10 ; 11; 12|] + in + print_endline @@ to_string Kind_2 k2 ; + print_endline @@ to_string Kind_3 k3 ; + print_endline @@ to_string Kind_4 k4 ; + print_endline @@ to_string Kind_6 k6 ; + print_endline @@ to_string Kind_12 k12 + + +let () = debug () +*) diff --git a/Utils/Zmap.ml b/Utils/Zmap.ml new file mode 100644 index 0000000..e3c5a85 --- /dev/null +++ b/Utils/Zmap.ml @@ -0,0 +1,4 @@ +(** Hash table where the keys are of type Zkey.t (tuples of integers) *) + +module Zmap = Hashtbl.Make(Zkey) +include Zmap