From ade45cd49d91eac3b92f94aac5a8809764559e04 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 13 Feb 2018 17:36:25 +0100 Subject: [PATCH] A lot of cleaning --- Basis/Contracted_shell.ml | 7 +- Basis/Contracted_shell_type.ml | 73 +++++++++++++++++++ Basis/Contracted_shell_type.mli | 24 ++++++ Basis/KinInt.ml | 10 ++- Basis/OneElectronRR.ml | 31 ++++---- Basis/Overlap.ml | 10 ++- Basis/TwoElectronRR.ml | 76 +++++++++---------- Basis/TwoElectronRRVectorized.ml | 82 +++++++++++---------- Nuclei/Nuclei.ml | 24 +++++- Nuclei/Xyz_ast.mli | 13 ++++ Nuclei/Xyz_parser.mly | 59 ++++++++++----- Utils/Angstrom.ml | 12 +++ Utils/Angstrom.mli | 10 +++ Utils/Bohr.ml | 12 +++ Utils/Bohr.mli | 10 +++ Utils/Constants.ml | 1 + Utils/Constants.mli | 8 ++ Utils/Coordinate.ml | 121 ++++++++++++++----------------- Utils/Coordinate.mli | 36 ++++----- Utils/Coordinate_type.ml | 26 +++++++ Utils/Coordinate_type.mli | 8 ++ Utils/Point.mli | 7 ++ 22 files changed, 450 insertions(+), 210 deletions(-) create mode 100644 Basis/Contracted_shell_type.ml create mode 100644 Basis/Contracted_shell_type.mli create mode 100644 Nuclei/Xyz_ast.mli create mode 100644 Utils/Angstrom.ml create mode 100644 Utils/Angstrom.mli create mode 100644 Utils/Bohr.ml create mode 100644 Utils/Bohr.mli create mode 100644 Utils/Constants.mli create mode 100644 Utils/Coordinate_type.ml create mode 100644 Utils/Coordinate_type.mli create mode 100644 Utils/Point.mli diff --git a/Basis/Contracted_shell.ml b/Basis/Contracted_shell.ml index a7f8fdc..5004d19 100644 --- a/Basis/Contracted_shell.ml +++ b/Basis/Contracted_shell.ml @@ -1,6 +1,7 @@ open Util open Constants open Contracted_shell_type +open Coordinate type t = Contracted_shell_type.t @@ -16,16 +17,14 @@ let with_index = Contracted_shell_type.with_index let powers a = a.powers let to_string s = - let coord = - Coordinate.to_Bohr s.center - in + let coord = s.center in let open Printf in (match s.totAngMom with | Angular_momentum.S -> sprintf "%3d " (s.index+1) | _ -> sprintf "%3d-%-3d" (s.index+1) (s.index+(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) ) ^ + (get X coord) (get Y coord) (get 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" " ") ) diff --git a/Basis/Contracted_shell_type.ml b/Basis/Contracted_shell_type.ml new file mode 100644 index 0000000..ab7d494 --- /dev/null +++ b/Basis/Contracted_shell_type.ml @@ -0,0 +1,73 @@ +open Util +open Constants + +type shell_contracted = { + expo : float array; + coef : float array; + center : Coordinate.t; + totAngMom : Angular_momentum.t; + size : int; + norm_coef : float array; + norm_coef_scale : float array; + index : int; + powers : Zkey.t array; +} + +type t = shell_contracted + +(** 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. + Returns, for each contracted function, an array of functions taking as + argument the [|x;y;z|] powers. +*) +let compute_norm_coef expo totAngMom = + let atot = + Angular_momentum.to_int totAngMom + in + let factor int_array = + let dfa = Array.map (fun j -> + ( float_of_int (1 lsl j) *. fact j) /. fact (j+j) + ) int_array + in + sqrt (dfa.(0) *.dfa.(1) *. dfa.(2)) + in + let expo = + if atot mod 2 = 0 then + Array.map (fun alpha -> + let alpha_2 = alpha +. alpha in + (alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2)) + ) expo + else + Array.map (fun alpha -> + let alpha_2 = alpha +. alpha in + (alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot) + ) expo + in + Array.map (fun x -> let f a = x *. (factor a) in f) expo + + +let make ~index ~expo ~coef ~center ~totAngMom = + assert (Array.length expo = Array.length coef); + assert (Array.length expo > 0); + let norm_coef_func = + compute_norm_coef expo totAngMom + in + let powers = + Angular_momentum.zkey_array (Angular_momentum.Singlet totAngMom) + in + let norm_coef = + Array.map (fun f -> f [| Angular_momentum.to_int totAngMom ; 0 ; 0 |]) norm_coef_func + in + let norm_coef_scale = + Array.map (fun a -> + (norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0) + ) powers + in + { index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ; + norm_coef_scale ; powers } + + +let with_index a i = { a with index = i } + diff --git a/Basis/Contracted_shell_type.mli b/Basis/Contracted_shell_type.mli new file mode 100644 index 0000000..1ec55fb --- /dev/null +++ b/Basis/Contracted_shell_type.mli @@ -0,0 +1,24 @@ +type shell_contracted = private { + expo : float array; + coef : float array; + center : Coordinate.t; + totAngMom : Angular_momentum.t; + size : int; + norm_coef : float array; + norm_coef_scale : float array; + index : int; + powers : Zkey.t array; +} + +type t = shell_contracted + + +val make : + index:int -> + expo:float array -> + coef:float array -> + center:Coordinate.t -> totAngMom:Angular_momentum.t -> t + + +(** Returns a copy of the contracted shell with a modified index *) +val with_index : t -> int -> t diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index 57106bc..526ea36 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -1,6 +1,7 @@ open Util open Constants open Lacaml.D +open Coordinate type t = Mat.t @@ -62,10 +63,15 @@ let contracted_class shell_a shell_b : float Zmap.t = [| a.(3) ; a.(4) ; a.(5) |] ) in let ov a b k = + let xyz = match k with + | 0 -> X + | 1 -> Y + | _ -> Z + in Overlap_primitives.hvrr (a, b) expo_inv - (Coordinate.coord center_ab k, - Coordinate.coord center_pa k) + (Coordinate.get xyz center_ab, + Coordinate.get xyz center_pa) in let f k = ov angMomA.(k) angMomB.(k) k diff --git a/Basis/OneElectronRR.ml b/Basis/OneElectronRR.ml index f517ba2..38f72e5 100644 --- a/Basis/OneElectronRR.ml +++ b/Basis/OneElectronRR.ml @@ -1,5 +1,6 @@ open Util open Constants +open Coordinate exception NullPair @@ -27,7 +28,7 @@ let hvrr_one_e (** Vertical recurrence relations *) let rec vrr angMom_a totAngMom_a = let ax,ay,az = angMom_a in - if (ax < 0) || (ay < 0) || (az < 0) then + if ax < 0 || ay < 0 || az < 0 then empty else match totAngMom_a with @@ -40,13 +41,13 @@ let hvrr_one_e let result = let am, amm, amxyz, xyz = match angMom_a with - | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, 0 - | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1 - | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2 + | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, X + | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, Y + | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, Z in if amxyz < 0 then empty else - let f1 = Coordinate.coord center_pa xyz - and f2 = expo_inv_p *. (Coordinate.coord center_pc xyz) + let f1 = Coordinate.get xyz center_pa + and f2 = expo_inv_p *. (Coordinate.get xyz center_pc) in if amxyz < 1 then let v1 = @@ -77,7 +78,7 @@ let hvrr_one_e and hrr angMom_a angMom_b totAngMom_a totAngMom_b = let bx,by,bz = angMom_b in - if (bx < 0) || (by < 0) || (bz < 0) then 0. + if bx < 0 || by < 0 || bz < 0 then 0. else match totAngMom_b with | 0 -> (vrr angMom_a totAngMom_a).(0) @@ -86,24 +87,24 @@ let hvrr_one_e and angMom_bx, angMom_by, angMom_bz = angMom_b in let bxyz, xyz = match angMom_b with - | (_,0,0) -> angMom_bx, 0 - | (_,_,0) -> angMom_by, 1 - | (_,_,_) -> angMom_bz, 2 + | (_,0,0) -> angMom_bx, X + | (_,_,0) -> angMom_by, Y + | (_,_,_) -> angMom_bz, Z in if (bxyz < 1) then 0. else let ap, bm = match xyz with - | 0 -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz) - | 1 -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz) - | _ -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1) + | X -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz) + | Y -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz) + | Z -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1) in let h1 = hrr ap bm (totAngMom_a+1) (totAngMom_b-1) in let f2 = - (Coordinate.coord center_ab xyz) + Coordinate.get xyz center_ab in - if (abs_float f2 < cutoff) then h1 else + if abs_float f2 < cutoff then h1 else let h2 = hrr angMom_a bm totAngMom_a (totAngMom_b-1) in diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index d65e22a..e40b4c3 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -1,6 +1,7 @@ open Util open Constants open Lacaml.D +open Coordinate type t = Mat.t @@ -56,10 +57,15 @@ let contracted_class shell_a shell_b : float Zmap.t = [| a.(3) ; a.(4) ; a.(5) |] ) in let f k = + let xyz = match k with + | 0 -> X + | 1 -> Y + | _ -> Z + in Overlap_primitives.hvrr (angMomA.(k), angMomB.(k)) expo_inv - (Coordinate.coord center_ab k, - Coordinate.coord center_pa k) + (Coordinate.get xyz center_ab, + Coordinate.get xyz center_pa) in let norm = norm_coef_scale.(i) in let integral = chop norm (fun () -> (f 0)*.(f 1)*.(f 2)) in diff --git a/Basis/TwoElectronRR.ml b/Basis/TwoElectronRR.ml index f4ed9da..9c85841 100644 --- a/Basis/TwoElectronRR.ml +++ b/Basis/TwoElectronRR.ml @@ -1,5 +1,6 @@ open Util open Constants +open Coordinate let debug=false @@ -45,9 +46,9 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let (x,y,z) = angMom_d in Printf.printf "%d %d %d\n" x y z; Printf.printf "%f %f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n" expo_b expo_d expo_inv_p expo_inv_q - (Coordinate.coord center_ab 0) (Coordinate.coord center_ab 1) (Coordinate.coord center_ab 2) - (Coordinate.coord center_cd 0) (Coordinate.coord center_cd 1) (Coordinate.coord center_cd 2) - (Coordinate.coord center_pq 0) (Coordinate.coord center_pq 1) (Coordinate.coord center_pq 2) + (get X center_ab) (get Y center_ab) (get Z center_ab) + (get X center_cd) (get Y center_cd) (get Z center_cd) + (get X center_pq) (get Y center_pq) (get Z center_pq) end; (** Vertical recurrence relations *) @@ -68,13 +69,13 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let result = let am, amm, amxyz, xyz = match angMom_a with - | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, 0 - | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1 - | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2 + | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, X + | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, Y + | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, Z in if amxyz < 0 then empty else - let f1 = expo_inv_p *. (Coordinate.coord center_pq xyz) - and f2 = expo_b *. expo_inv_p *. (Coordinate.coord center_ab xyz) + let f1 = expo_inv_p *. (Coordinate.get xyz center_pq) + and f2 = expo_b *. expo_inv_p *. (Coordinate.get xyz center_ab) in let result = Array.create_float maxsze in if amxyz < 1 then @@ -140,24 +141,24 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) (angMom_ax-1, angMom_ay, angMom_az), (angMom_cx-1, angMom_cy, angMom_cz), (angMom_cx-2, angMom_cy, angMom_cz), - angMom_ax,angMom_cx-1, 0 + angMom_ax,angMom_cx-1, X | (_,_,0) -> (angMom_ax, angMom_ay-1, angMom_az), (angMom_cx, angMom_cy-1, angMom_cz), (angMom_cx, angMom_cy-2, angMom_cz), - angMom_ay,angMom_cy-1, 1 + angMom_ay,angMom_cy-1, Y | _ -> (angMom_ax, angMom_ay, angMom_az-1), (angMom_cx, angMom_cy, angMom_cz-1), (angMom_cx, angMom_cy, angMom_cz-2), - angMom_az,angMom_cz-1, 2 + angMom_az,angMom_cz-1, Z in if cmxyz < 0 then empty else let f1 = - -. expo_d *. expo_inv_q *. (Coordinate.coord center_cd xyz) + -. expo_d *. expo_inv_q *. (Coordinate.get xyz center_cd) and f2 = - expo_inv_q *. (Coordinate.coord center_pq xyz) + expo_inv_q *. (Coordinate.get xyz center_pq) in let result = Array.make maxsze 0. in if ( (abs_float f1 > cutoff) || (abs_float f2 > cutoff) ) then @@ -213,6 +214,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and hrr0 angMom_a angMom_b angMom_c totAngMom_a totAngMom_b totAngMom_c = + (* if debug then begin let angMom_ax, angMom_ay, angMom_az = angMom_a @@ -224,6 +226,7 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) angMom_bx angMom_by angMom_bz angMom_cx angMom_cy angMom_cz end; + *) match totAngMom_b with | 0 -> (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0) @@ -231,15 +234,15 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) let angMom_ax, angMom_ay, angMom_az = angMom_a in let ap, xyz = match angMom_b with - | (1,_,_) -> (angMom_ax+1,angMom_ay,angMom_az), 0 - | (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), 1 - | _ -> (angMom_ax,angMom_ay,angMom_az+1), 2 + | (1,_,_) -> (angMom_ax+1,angMom_ay,angMom_az), X + | (_,1,_) -> (angMom_ax,angMom_ay+1,angMom_az), Y + | _ -> (angMom_ax,angMom_ay,angMom_az+1), Z in let v1 = vrr ap angMom_c (totAngMom_a+1) totAngMom_c in let f2 = - (Coordinate.coord center_ab xyz) + (Coordinate.get xyz center_ab) in if (abs_float f2 < cutoff) then v1.(0) else let v2 = @@ -251,22 +254,22 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and angMom_bx, angMom_by, angMom_bz = angMom_b in let bxyz, xyz = match angMom_b with - | (_,0,0) -> angMom_bx, 0 - | (_,_,0) -> angMom_by, 1 - | (_,_,_) -> angMom_bz, 2 + | (_,0,0) -> angMom_bx, X + | (_,_,0) -> angMom_by, Y + | (_,_,_) -> angMom_bz, Z in if (bxyz < 1) then 0. else let ap, bm = match xyz with - | 0 -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz) - | 1 -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz) - | _ -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1) + | X -> (angMom_ax+1,angMom_ay,angMom_az),(angMom_bx-1,angMom_by,angMom_bz) + | Y -> (angMom_ax,angMom_ay+1,angMom_az),(angMom_bx,angMom_by-1,angMom_bz) + | Z -> (angMom_ax,angMom_ay,angMom_az+1),(angMom_bx,angMom_by,angMom_bz-1) in let h1 = hrr0 ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c in let f2 = - (Coordinate.coord center_ab xyz) + (Coordinate.get xyz center_ab) in if (abs_float f2 < cutoff) then h1 else let h2 = @@ -274,23 +277,10 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) in h1 +. f2 *. h2 + and hrr angMom_a angMom_b angMom_c angMom_d totAngMom_a totAngMom_b totAngMom_c totAngMom_d = - if debug then - begin - let angMom_ax, angMom_ay, angMom_az = angMom_a in - let angMom_bx, angMom_by, angMom_bz = angMom_b in - let angMom_cx, angMom_cy, angMom_cz = angMom_c in - let angMom_dx, angMom_dy, angMom_dz = angMom_d in - Printf.printf "hrr : %d %d %d %d : %d %d %d %d %d %d %d %d %d %d %d %d\n" - totAngMom_a totAngMom_b totAngMom_c totAngMom_d - angMom_ax angMom_ay angMom_az - angMom_bx angMom_by angMom_bz - angMom_cx angMom_cy angMom_cz - angMom_dx angMom_dy angMom_dz - end; - match (totAngMom_b, totAngMom_d) with | (_,0) -> if (totAngMom_b = 0) then (vrr angMom_a angMom_c totAngMom_a totAngMom_c).(0) @@ -301,14 +291,14 @@ let rec hvrr_two_e (angMom_a, angMom_b, angMom_c, angMom_d) and (angMom_dx, angMom_dy, angMom_dz) = angMom_d in let cp, dm, xyz = match angMom_d with - | (_,0,0) -> (angMom_cx+1, angMom_cy, angMom_cz), (angMom_dx-1, angMom_dy, angMom_dz), 0 - | (_,_,0) -> (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), 1 - | _ -> (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), 2 + | (_,0,0) -> (angMom_cx+1, angMom_cy, angMom_cz), (angMom_dx-1, angMom_dy, angMom_dz), X + | (_,_,0) -> (angMom_cx, angMom_cy+1, angMom_cz), (angMom_dx, angMom_dy-1, angMom_dz), Y + | _ -> (angMom_cx, angMom_cy, angMom_cz+1), (angMom_dx, angMom_dy, angMom_dz-1), Z in let h1 = hrr angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) in - let f2 = Coordinate.coord center_cd xyz in + let f2 = Coordinate.get xyz center_cd in if (abs_float f2 < cutoff) then h1 else let h2 = hrr angMom_a angMom_b angMom_c dm totAngMom_a totAngMom_b totAngMom_c (totAngMom_d-1) @@ -372,7 +362,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q raise NullQuartet; let center_pq = - Coordinate.(sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center) + sp.(ab).ShellPair.center |- sq.(cd).ShellPair.center in let norm_pq_sq = Coordinate.dot center_pq center_pq diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index 6b046d4..5efaafa 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -1,6 +1,7 @@ open Util open Lacaml.D open Bigarray +open Coordinate let cutoff = Constants.cutoff let cutoff2 = cutoff *. cutoff @@ -42,15 +43,15 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let result = let am, amm, amxyz, xyz = match angMom_a with - | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, 0 - | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, 1 - | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, 2 + | (x,0,0) -> (x-1,0,0),(x-2,0,0), x-1, X + | (x,y,0) -> (x,y-1,0),(x,y-2,0), y-1, Y + | (x,y,z) -> (x,y,z-1),(x,y,z-2), z-1, Z in if amxyz < 0 then None else begin - let cab = Coordinate.coord center_ab xyz in + let cab = Coordinate.get xyz center_ab in let v1_top, p1_top = if abs_float cab < cutoff then None, @@ -89,7 +90,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) for l=0 to np-1 do for k=0 to nq-1 do result.(l).(k) <- result.(l).(k) - +. expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k) + +. expo_inv_p.(l) *. (center_pq xyz).(l).(k) *. p0.(l).(k) done done end @@ -126,7 +127,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) for k=0 to nq-1 do result.(l).(k) <- result.(l).(k) +. - expo_inv_p.(l) *. center_pq.(xyz).(l).(k) *. p0.(l).(k) +. + expo_inv_p.(l) *. (center_pq xyz).(l).(k) *. p0.(l).(k) +. f *. (v1.(l).(k) +. v2.(l).(k) *. expo_inv_p.(l)) done done @@ -157,24 +158,24 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) (aax, aay, aaz-1), (acx, acy, acz-1), (acx, acy, acz-2), - aaz, acz, 2 + aaz, acz, Z else if (acy > 0) then (aax, aay-1,aaz), (acx, acy-1,acz), (acx, acy-2,acz), - aay,acy, 1 + aay,acy, Y else (aax-1,aay,aaz), (acx-1,acy,acz), (acx-2,acy,acz), - aax,acx, 0 + aax,acx, X in (* let result = Array.make_matrix np nq 0. in *) let do_compute = ref false in let v1 = - let f = -. (Coordinate.coord center_cd xyz) in + let f = -. (Coordinate.get xyz center_cd) in let f1 = Array.init nq (fun k -> @@ -203,7 +204,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let f2 = Array.init np (fun l -> Array.init nq (fun k -> - let x = expo_inv_q.(k) *. center_pq.(xyz).(l).(k) in + let x = expo_inv_q.(k) *. (center_pq xyz).(l).(k) in if (!do_compute) then x else (if abs_float x > cutoff then do_compute := true ; x) ) ) @@ -382,11 +383,11 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) let (aax, aay, aaz) = angMom_a in let ap, xyz = match angMom_b with - | (_,_,1) -> (aax,aay,aaz+1), 2 - | (_,1,_) -> (aax,aay+1,aaz), 1 - | (_,_,_) -> (aax+1,aay,aaz), 0 + | (_,_,1) -> (aax,aay,aaz+1), Z + | (_,1,_) -> (aax,aay+1,aaz), Y + | (_,_,_) -> (aax+1,aay,aaz), X in - let f = Coordinate.coord center_ab xyz in + let f = Coordinate.get xyz center_ab in let v1 = match vrr_v 0 ap angMom_c (totAngMom_a+1) totAngMom_c with | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix @@ -404,23 +405,23 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) and (abx, aby, abz) = angMom_b in let bxyz, xyz = match angMom_b with - | (0,0,_) -> abz, 2 - | (0,_,_) -> aby, 1 - | _ -> abx, 0 + | (0,0,_) -> abz, Z + | (0,_,_) -> aby, Y + | _ -> abx, X in if (bxyz < 1) then 0. else let ap, bm = match xyz with - | 0 -> (aax+1,aay,aaz),(abx-1,aby,abz) - | 1 -> (aax,aay+1,aaz),(abx,aby-1,abz) - | _ -> (aax,aay,aaz+1),(abx,aby,abz-1) + | X -> (aax+1,aay,aaz),(abx-1,aby,abz) + | Y -> (aax,aay+1,aaz),(abx,aby-1,abz) + | Z -> (aax,aay,aaz+1),(abx,aby,abz-1) in let h1 = hrr0_v ap bm angMom_c (totAngMom_a+1) (totAngMom_b-1) totAngMom_c in - let f = (Coordinate.coord center_ab xyz) in - if (abs_float f < cutoff) then h1 else + let f = Coordinate.get xyz center_ab in + if abs_float f < cutoff then h1 else let h2 = hrr0_v angMom_a bm angMom_c totAngMom_a (totAngMom_b-1) totAngMom_c in @@ -430,7 +431,7 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) totAngMom_a totAngMom_b totAngMom_c totAngMom_d = match (totAngMom_b, totAngMom_d) with - | (_,0) -> if (totAngMom_b = 0) then + | (_,0) -> if totAngMom_b = 0 then begin match vrr_v 0 angMom_a angMom_c totAngMom_a totAngMom_c with | Some matrix -> Array.fold_left (fun accu c -> accu +. Array.fold_left (+.) 0. c) 0. matrix @@ -443,15 +444,15 @@ let hvrr_two_e_vector (angMom_a, angMom_b, angMom_c, angMom_d) and (adx, ady, adz) = angMom_d in let cp, dm, xyz = match angMom_d with - | (_,0,0) -> (acx+1, acy, acz), (adx-1, ady, adz), 0 - | (_,_,0) -> (acx, acy+1, acz), (adx, ady-1, adz), 1 - | _ -> (acx, acy, acz+1), (adx, ady, adz-1), 2 + | (_,0,0) -> (acx+1, acy, acz), (adx-1, ady, adz), X + | (_,_,0) -> (acx, acy+1, acz), (adx, ady-1, adz), Y + | _ -> (acx, acy, acz+1), (adx, ady, adz-1), Z in let h1 = hrr_v angMom_a angMom_b cp dm totAngMom_a totAngMom_b (totAngMom_c+1) (totAngMom_d-1) in - let f = (Coordinate.coord center_cd xyz) in - if (abs_float f < cutoff) then + let f = Coordinate.get xyz center_cd in + if abs_float f < cutoff then h1 else let h2 = @@ -570,7 +571,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in let center_pq = - Coordinate.(sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center) + sp.(i-1).ShellPair.center |- sq.(j-1).ShellPair.center in let norm_pq_sq = Coordinate.dot center_pq center_pq @@ -608,6 +609,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let norm_coef_scale_p = shell_p.ContractedShellPair.norm_coef_scale in let center_pq = + let result = Array.init 3 (fun xyz -> Array.init np (fun ab -> let shell_ab = sp.(ab) in @@ -615,16 +617,20 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q let shell_cd = sq.(cd) in let cpq = - Coordinate.(shell_ab.ShellPair.center |- shell_cd.ShellPair.center) + shell_ab.ShellPair.center |- shell_cd.ShellPair.center in match xyz with - | 0 -> Coordinate.x cpq; - | 1 -> Coordinate.y cpq; - | 2 -> Coordinate.z cpq; + | 0 -> Coordinate.get X cpq; + | 1 -> Coordinate.get Y cpq; + | 2 -> Coordinate.get Z cpq; | _ -> assert false ) ) ) + in function + | X -> result.(0) + | Y -> result.(1) + | Z -> result.(2) in let zero_m_array = let result = @@ -642,9 +648,11 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q expo_inv_p.(ab) +. expo_inv_q.(cd) in let norm_pq_sq = - center_pq.(0).(ab).(cd) *. center_pq.(0).(ab).(cd) +. - center_pq.(1).(ab).(cd) *. center_pq.(1).(ab).(cd) +. - center_pq.(2).(ab).(cd) *. center_pq.(2).(ab).(cd) + let x = (center_pq X).(ab).(cd) + and y = (center_pq Y).(ab).(cd) + and z = (center_pq Z).(ab).(cd) + in + x *. x +. y *. y +. z *. z in zero_m ~maxm ~expo_pq_inv ~norm_pq_sq diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index c1babb6..cab3cad 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -1,3 +1,5 @@ +open Xyz_ast + type t = (Element.t * Coordinate.t) array let of_xyz_file ~filename = @@ -5,7 +7,21 @@ let of_xyz_file ~filename = let ic = open_in filename in Lexing.from_channel ic in - Xyz_parser.input Nuclei_lexer.read_all lexbuf + let data = + Xyz_parser.input Nuclei_lexer.read_all lexbuf + in + + let len = List.length data.nuclei in + if len <> data.number_of_atoms then + Printf.sprintf "Error: expected %d atoms but %d read" + data.number_of_atoms len + |> failwith; + + List.map (fun nucleus -> + nucleus.element, Coordinate.angstrom_to_bohr nucleus.coord + ) data.nuclei + |> Array.of_list + let of_zmt_file ~filename = @@ -21,7 +37,7 @@ let of_zmt_file ~filename = in aux [] |> Zmatrix.of_string |> Zmatrix.to_xyz - |> Array.map (fun (e,x,y,z) -> (e, Coordinate.of_3_floats x y z )) + |> Array.map (fun (e,x,y,z) -> (e, (Angstrom.make {Point.x ; y ; z} ))) let to_string atoms = @@ -36,11 +52,11 @@ let to_string atoms = " ^ (Array.mapi (fun i (e, coord) -> let coord = - Coordinate.to_Angstrom coord + Coordinate.bohr_to_angstrom coord in Printf.sprintf " %5d %5d %5s %12.6f %12.6f %12.6f" (i+1) (Element.to_int e) (Element.to_string e) - (Coordinate.x coord) (Coordinate.y coord) (Coordinate.z coord) + coord.Angstrom.x coord.Angstrom.y coord.Angstrom.z ) atoms |> Array.to_list |> String.concat "\n" ) ^ diff --git a/Nuclei/Xyz_ast.mli b/Nuclei/Xyz_ast.mli new file mode 100644 index 0000000..c1fc73f --- /dev/null +++ b/Nuclei/Xyz_ast.mli @@ -0,0 +1,13 @@ +type nucleus = + { + element: Element.t ; + coord: Angstrom.t; + } + +type xyz_file = + { + number_of_atoms : int; + file_title : string; + nuclei : nucleus list; + } + diff --git a/Nuclei/Xyz_parser.mly b/Nuclei/Xyz_parser.mly index 39bdf2e..e9c01d1 100644 --- a/Nuclei/Xyz_parser.mly +++ b/Nuclei/Xyz_parser.mly @@ -2,6 +2,24 @@ %{ exception InputError of string +let make_angstrom x y z = + Angstrom.make { + Point. + x ; y ; z + } + +let output_of f x y z = + let a = make_angstrom x y z in + fun e -> + { + Xyz_ast. + element = f e; + coord = a ; + } + +let output_of_string = output_of Element.of_string +let output_of_int = output_of Element.of_int + %} %token EOL @@ -12,54 +30,59 @@ exception InputError of string %token EOF %start input -%type <(Element.t * Coordinate.t) array> input +%type input %% /* Grammar rules and actions follow */ input: | integer title atoms_xyz { - let len = List.length $3 in - if len <> $1 then - let error_msg = Printf.sprintf "%d atoms entered, expected %d" len $1 in - raise (InputError error_msg) - else - Array.of_list $3 - } + { + number_of_atoms = $1; + file_title = $2; + nuclei = $3; + } + } +; + integer: | INTEGER EOL { $1 } | INTEGER SPACE EOL { $1 } | SPACE INTEGER EOL { $2 } | SPACE INTEGER SPACE EOL { $2 } +; title: | title_list EOL { $1 } +; text: | WORD { $1 } | SPACE { $1 } | FLOAT { (string_of_float $1)} | INTEGER { (string_of_int $1)} +; title_list: | { "" } | title_list text { ($1 ^ $2) } +; atoms_xyz: | atoms_list EOL { List.rev $1 } | atoms_list EOF { List.rev $1 } - +; atoms_list: | { [] } - | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_string $2, Coordinate.of_3_floats $4 $6 $8 `Angstrom ) :: $1 } - | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_string $2, Coordinate.of_3_floats $4 $6 $8 `Angstrom ) :: $1 } - | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_int $2, Coordinate.of_3_floats $4 $6 $8 `Angstrom ) :: $1 } - | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_int $2, Coordinate.of_3_floats $4 $6 $8 `Angstrom ) :: $1 } - | atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_string $3, Coordinate.of_3_floats $5 $7 $9 `Angstrom ) :: $1 } - | atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_string $3, Coordinate.of_3_floats $5 $7 $9 `Angstrom ) :: $1 } - | atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { (Element.of_int $3, Coordinate.of_3_floats $5 $7 $9 `Angstrom ) :: $1 } - | atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { (Element.of_int $3, Coordinate.of_3_floats $5 $7 $9 `Angstrom ) :: $1 } - + | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_string $4 $6 $8 $2 :: $1 } + | atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_string $4 $6 $8 $2 :: $1 } + | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_int $4 $6 $8 $2 :: $1 } + | atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_int $4 $6 $8 $2 :: $1 } + | atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_string $5 $7 $9 $3 :: $1 } + | atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_string $5 $7 $9 $3 :: $1 } + | atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_int $5 $7 $9 $3 :: $1 } + | atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_int $5 $7 $9 $3 :: $1 } +; diff --git a/Utils/Angstrom.ml b/Utils/Angstrom.ml new file mode 100644 index 0000000..b553427 --- /dev/null +++ b/Utils/Angstrom.ml @@ -0,0 +1,12 @@ +type t = { + x : float ; + y : float ; + z : float ; +} + + +let make { Point.x ; y ; z } = { x ; y ; z } +(* = %identity *) + + + diff --git a/Utils/Angstrom.mli b/Utils/Angstrom.mli new file mode 100644 index 0000000..77aa9e2 --- /dev/null +++ b/Utils/Angstrom.mli @@ -0,0 +1,10 @@ +type t = private { + x : float ; + y : float ; + z : float ; +} + + +val make : Point.t -> t + + diff --git a/Utils/Bohr.ml b/Utils/Bohr.ml new file mode 100644 index 0000000..b553427 --- /dev/null +++ b/Utils/Bohr.ml @@ -0,0 +1,12 @@ +type t = { + x : float ; + y : float ; + z : float ; +} + + +let make { Point.x ; y ; z } = { x ; y ; z } +(* = %identity *) + + + diff --git a/Utils/Bohr.mli b/Utils/Bohr.mli new file mode 100644 index 0000000..77aa9e2 --- /dev/null +++ b/Utils/Bohr.mli @@ -0,0 +1,10 @@ +type t = private { + x : float ; + y : float ; + z : float ; +} + + +val make : Point.t -> t + + diff --git a/Utils/Constants.ml b/Utils/Constants.ml index 1f855f7..9d8f9f9 100644 --- a/Utils/Constants.ml +++ b/Utils/Constants.ml @@ -8,3 +8,4 @@ let pi_inv = 1. /. pi let two_over_sq_pi = 2. /. (sqrt pi) let a0 = 0.529_177_210_67 +let a0_inv = 1. /. a0 diff --git a/Utils/Constants.mli b/Utils/Constants.mli new file mode 100644 index 0000000..48541b6 --- /dev/null +++ b/Utils/Constants.mli @@ -0,0 +1,8 @@ +val cutoff : float +val pi : float +val sq_pi : float +val sq_pi_over_two : float +val pi_inv : float +val two_over_sq_pi : float +val a0 : float +val a0_inv : float diff --git a/Utils/Coordinate.ml b/Utils/Coordinate.ml index eaabfa5..9daa642 100644 --- a/Utils/Coordinate.ml +++ b/Utils/Coordinate.ml @@ -1,89 +1,76 @@ -open Util +type bohr = Bohr.t +type angstrom = Angstrom.t -type t = - | Bohr of (float * float * float) - | Angstrom of (float * float * float) +type t = bohr -let zero = Bohr (0., 0., 0.) +type axis = X | Y | Z -let of_float_triplet (x,y,z) = function -| `Bohr -> Bohr (x,y,z) -| `Angstrom -> Angstrom (x,y,z) +let a_to_b a = Constants.a0_inv *. a +let b_to_a b = Constants.a0 *. b -let of_3_floats x y z = - of_float_triplet (x,y,z) +let bohr_to_angstrom { Bohr.x ; y ; z } = + Angstrom.make + { + Point. + x = b_to_a x ; + y = b_to_a y ; + z = b_to_a z ; + } -let to_string t = - let result (x,y,z) = - (string_of_float x)^" "^(string_of_float y)^" "^(string_of_float z) - in - match t with - | Bohr x -> (result x) ^ " Bohr" - | Angstrom x -> (result x) ^ " Angstrom" +let angstrom_to_bohr { Angstrom.x ; y ; z } = + Bohr.make + { + Point. + x = a_to_b x ; + y = a_to_b y ; + z = a_to_b z ; + } -let extract_float_tuple = function - | Bohr a - | Angstrom a -> a +let zero = + Bohr.make { Point.x = 0. ; y = 0. ; z = 0. } (** Linear algebra *) -let (|.) s a = - match a with - | Bohr (x,y,z) -> Bohr ( s*.x, s*.y, s*.z ) - | Angstrom (x,y,z) -> Angstrom ( s*.x, s*.y, s*.z ) +let ( |. ) s { Bohr.x ; y ; z } = + Bohr.make { + Point. + x = s *. x ; + y = s *. y ; + z = s *. z ; + } -let to_Angstrom = function - | Angstrom a -> Angstrom a - | Bohr a -> Angstrom (Constants.a0 |. Bohr a |> extract_float_tuple) +let ( |+ ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } = + Bohr.make { + Point. + x = x1 +. x2 ; + y = y1 +. y2 ; + z = z1 +. z2 ; + } -let to_Bohr = function - | Angstrom a -> Bohr (1./.Constants.a0 |. Angstrom a |> extract_float_tuple) - | Bohr a -> Bohr a -let (|-), (|+) = - let rec op f p q = - match (p, q) with - | (Angstrom a, Angstrom b) -> Angstrom (f a b) - | (Bohr a, Bohr b) -> Bohr (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 (x,y,z) (x',y',z') -> ( x-.x', y-.y', z-.z' )) , - op (fun (x,y,z) (x',y',z') -> ( x+.x', y+.y', z+.z' )) - ) +let ( |- ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } = + Bohr.make { + Point. + x = x1 -. x2 ; + y = y1 -. y2 ; + z = z1 -. z2 ; + } let neg a = -1. |. a -let rec dot p q = - match (p,q) with - | Bohr (x,y,z), Bohr (x',y',z') -> x*.x' +. y*.y' +. z*.z' - | _ -> dot (to_Bohr p) (to_Bohr q) + +let dot { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } = + x1 *. x2 +. y1 *. y2 +. z1 *. z2 let norm u = - sqrt @@ dot u u + sqrt ( dot u u ) -let rec to_tuple a = - to_Bohr a |> extract_float_tuple - -let x a = - let (result, _, _) = extract_float_tuple @@ to_Bohr a in - result - -let y a = - let (_, result, _) = extract_float_tuple @@ to_Bohr a in - result - -let z a = - let (_, _, result) = extract_float_tuple @@ to_Bohr a in - result - -let coord a = function - | 0 -> x a - | 1 -> y a - | 2 -> z a - | _ -> raise (Invalid_argument "Coordinate") - +let get axis { Bohr.x ; y ; z } = + match axis with + | X -> x + | Y -> y + | Z -> z diff --git a/Utils/Coordinate.mli b/Utils/Coordinate.mli index 13df088..dee7a71 100644 --- a/Utils/Coordinate.mli +++ b/Utils/Coordinate.mli @@ -1,19 +1,19 @@ -type t -val to_Angstrom : t -> t -val to_Bohr : t -> t -val zero : t -val of_float_triplet : (float * float * float) -> [< `Angstrom | `Bohr ] -> t -val of_3_floats : float -> float -> float -> [< `Angstrom | `Bohr ] -> t -val ( |. ) : float -> t -> t -val ( |- ) : t -> t -> t -val ( |+ ) : t -> t -> t -val neg : t -> t -val dot : t -> t -> float -val norm : t -> float -val to_string : t -> string -val to_tuple : t -> (float * float * float) -val x : t -> float -val y : t -> float -val z : t -> float -val coord : t -> int -> float +type bohr = Bohr.t +type angstrom = Angstrom.t + +type t = bohr +type axis = X | Y | Z + +val bohr_to_angstrom : bohr -> angstrom +val angstrom_to_bohr : angstrom -> bohr + +val zero : bohr + +val ( |. ) : float -> bohr -> bohr +val ( |+ ) : bohr -> bohr -> bohr +val ( |- ) : bohr -> bohr -> bohr +val neg : bohr -> bohr +val dot : bohr -> bohr -> float +val norm : bohr -> float +val get : axis -> bohr -> float diff --git a/Utils/Coordinate_type.ml b/Utils/Coordinate_type.ml new file mode 100644 index 0000000..b2af9a8 --- /dev/null +++ b/Utils/Coordinate_type.ml @@ -0,0 +1,26 @@ +type point = Point.t + +type bohr = Bohr.t +type angstrom = Angstrom.t + +let a_to_b a = Constants.a0_inv *. a +let b_to_a b = Constants.a0 *. b + +let bohr_to_angstrom { Bohr.x ; y ; z } = + Angstrom.make + { + Point. + x = b_to_a x ; + y = b_to_a y ; + z = b_to_a z ; + } + +let angstrom_to_bohr { Angstrom.x ; y ; z } = + Bohr.make + { + Point. + x = a_to_b x ; + y = a_to_b y ; + z = a_to_b z ; + } + diff --git a/Utils/Coordinate_type.mli b/Utils/Coordinate_type.mli new file mode 100644 index 0000000..e3f1c00 --- /dev/null +++ b/Utils/Coordinate_type.mli @@ -0,0 +1,8 @@ +type bohr = Bohr.t +type angstrom = Angstrom.t + +val bohr_to_angstrom : bohr -> angstrom +val angstrom_to_bohr : angstrom -> bohr + + + diff --git a/Utils/Point.mli b/Utils/Point.mli new file mode 100644 index 0000000..78237ef --- /dev/null +++ b/Utils/Point.mli @@ -0,0 +1,7 @@ +type t = { + x : float ; + y : float ; + z : float ; +} + +