From b5b3e1fc98be03efb968a2d0456d377c77c64e1c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 26 Mar 2018 19:15:09 +0200 Subject: [PATCH] Added ERI type --- Basis/ERI.ml | 107 +++++++++++++++++++++++++++++++++++--------- Basis/ERI.mli | 1 - HartreeFock/Fock.ml | 4 +- 3 files changed, 89 insertions(+), 23 deletions(-) diff --git a/Basis/ERI.ml b/Basis/ERI.ml index c517f84..0c11e12 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -4,16 +4,88 @@ open Util open Constants open Bigarray -type t = (float, float32_elt, fortran_layout) Bigarray.Genarray.t +let max_ao = 1 lsl 14 -let get ~r1 ~r2 t = - let i,k = r1 - and j,l = r2 - in - t.{i,k,j,l} +type index_pair = { first : int ; second : int } -let get_chem t i j k l = get ~r1:(i,j) ~r2:(k,l) t -let get_phys t i j k l = get ~r1:(i,k) ~r2:(j,l) t +type t = + | Dense of (float, float64_elt, fortran_layout) Bigarray.Genarray.t + | Sparse of (int, float) Hashtbl.t + + +let key_of_indices ~r1 ~r2 = + let { first=i ; second=k } = r1 and { first=j ; second=l } = r2 in + let i,k = if i<=k then i,k else k,i + and j,l = if j<=l then j,l else l,j in + let i,k,j,l = if k<=l then i,k,j,l else j,l,i,k in + ((((((i lsl 15) lor k) lsl 15) lor j) lsl 15) lor l) + + +let get ~r1 ~r2 = function + | Dense t -> let { first=i ; second=k } = r1 and { first=j ; second=l } = r2 in t.{i,j,k,l} + | Sparse t -> let key = key_of_indices ~r1 ~r2 in + try Hashtbl.find t key + with Not_found -> 0. + + +let set ~r1 ~r2 ~value = function + | Dense t -> let { first=i ; second=k } = r1 and { first=j ; second=l } = r2 in + t.{i,j,k,l} <- value; + t.{k,j,i,l} <- value; + t.{i,l,k,j} <- value; + t.{k,l,i,j} <- value; + t.{j,i,l,k} <- value; + t.{j,k,l,i} <- value; + t.{l,i,j,k} <- value; + t.{l,k,j,i} <- value; + + | Sparse t -> let key = key_of_indices ~r1 ~r2 in + Hashtbl.replace t key value + + +let increment ~r1 ~r2 ~value = function + | Dense t -> let { first=i ; second=k } = r1 and { first=j ; second=l } = r2 in + t.{i,j,k,l} <- t.{i,j,k,l} +. value; + t.{k,j,i,l} <- t.{k,j,i,l} +. value; + t.{i,l,k,j} <- t.{i,l,k,j} +. value; + t.{k,l,i,j} <- t.{k,l,i,j} +. value; + t.{j,i,l,k} <- t.{j,i,l,k} +. value; + t.{j,k,l,i} <- t.{j,k,l,i} +. value; + t.{l,i,j,k} <- t.{l,i,j,k} +. value; + t.{l,k,j,i} <- t.{l,k,j,i} +. value; + | Sparse t -> let key = key_of_indices ~r1 ~r2 in + let old_value = + try Hashtbl.find t key + with Not_found -> 0. + in + Hashtbl.replace t key (old_value +. value) + + +let create = function + | `Dense n -> + let eri_array = + Genarray.create Float64 fortran_layout [| n ; n ; n ; n|] + in + Genarray.fill eri_array 0.; + Dense eri_array + + | `Sparse n -> + let eri_array = Hashtbl.create (n*n+13) in + Hashtbl.add eri_array (-1) (float_of_int n); + Sparse eri_array + +let size = function + | Dense t -> Genarray.nth_dim t 3 + | Sparse t -> Hashtbl.find t (-1) |> int_of_float + +(** TODO : remove epsilons *) + + +let get_chem t i j k l = get ~r1:{ first=i ; second=j } ~r2:{ first=k ; second=l } t +let get_phys t i j k l = get ~r1:{ first=i ; second=k } ~r2:{ first=j ; second=l } t + +let set_chem t i j k l value = set ~r1:{ first=i ; second=j } ~r2:{ first=k ; second=l } ~value t +let set_phys t i j k l value = set ~r1:{ first=i ; second=k } ~r2:{ first=j ; second=l } ~value t module Am = AngularMomentum @@ -140,9 +212,11 @@ let of_basis basis = (* 4D data initialization *) let eri_array = - Genarray.create Float32 fortran_layout [| n ; n ; n ; n|] + (* + create (`Dense n) + *) + create (`Sparse n) in - Genarray.fill eri_array 0.; (* Compute ERIs *) @@ -229,14 +303,7 @@ let of_basis basis = let value = Zmap.find cls key in - eri_array.{i_c,k_c,j_c,l_c} <- value; - eri_array.{j_c,k_c,i_c,l_c} <- value; - eri_array.{i_c,l_c,j_c,k_c} <- value; - eri_array.{j_c,l_c,i_c,k_c} <- value; - eri_array.{k_c,i_c,l_c,j_c} <- value; - eri_array.{k_c,j_c,l_c,i_c} <- value; - eri_array.{l_c,i_c,k_c,j_c} <- value; - eri_array.{l_c,j_c,k_c,i_c} <- value; + set_chem eri_array i_c j_c k_c l_c value; if (abs_float value > cutoff) then (inn := !inn + 1; ) @@ -259,11 +326,11 @@ let of_basis basis = let to_file ~filename eri_array = let oc = open_out filename in (* Print ERIs *) - for l_c=1 to (Genarray.nth_dim eri_array 3) do + for l_c=1 to size eri_array do for k_c=1 to l_c do for j_c=1 to l_c do for i_c=1 to k_c do - let value = eri_array.{i_c,j_c,k_c,l_c} in + let value = get_phys eri_array i_c j_c k_c l_c in if (abs_float value > cutoff) then Printf.fprintf oc " %5d %5d %5d %5d%20.15f\n" i_c j_c k_c l_c value; done; diff --git a/Basis/ERI.mli b/Basis/ERI.mli index 640e9d5..c906c63 100644 --- a/Basis/ERI.mli +++ b/Basis/ERI.mli @@ -1,5 +1,4 @@ type t -val get : r1:int * int -> r2:int * int -> t -> float val get_chem : t -> int -> int -> int -> int -> float val get_phys : t -> int -> int -> int -> int -> float diff --git a/HartreeFock/Fock.ml b/HartreeFock/Fock.ml index 62646eb..bb015b6 100644 --- a/HartreeFock/Fock.ml +++ b/HartreeFock/Fock.ml @@ -23,8 +23,8 @@ let make ~density simulation = if abs_float p > epsilon then for mu = 1 to nu do m_F.{mu,nu} <- m_F.{mu,nu} +. p *. - (ERI.get_chem m_G mu lambda nu sigma - -. 0.5 *. ERI.get_chem m_G mu lambda sigma nu) + (ERI.get_phys m_G mu lambda nu sigma + -. 0.5 *. ERI.get_phys m_G mu lambda sigma nu) done done done