diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 44337fa..d2828eb 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -3,7 +3,7 @@ open Util open Constants -type t = FourIdxStorage.t +include FourIdxStorage module Am = AngularMomentum @@ -16,14 +16,6 @@ module Cspc = ContractedShellPairCouple module Fis = FourIdxStorage -let get_chem = Fis.get_chem -let get_phys = Fis.get_phys - -let set_chem = Fis.set_chem -let set_phys = Fis.set_phys - -let to_file = Fis.to_file - let cutoff = integrals_cutoff diff --git a/Basis/ERI.mli b/Basis/ERI.mli index fded960..67206cf 100644 --- a/Basis/ERI.mli +++ b/Basis/ERI.mli @@ -1,6 +1,6 @@ (** Electron repulsion intergals. *) -type t +include module type of FourIdxStorage (* val filter_atomic_shell_pairs : ?cutoff:float -> AtomicShellPair.t list -> AtomicShellPair.t list @@ -13,16 +13,7 @@ val filter_contracted_shell_pairs : ?cutoff:float -> ContractedShellPair.t list val class_of_contracted_shell_pair_couple : ContractedShellPairCouple.t -> float Zmap.t (** Computes all the ERI of the class built from a couple of contracted shell pairs. *) - -val get_chem : t -> int -> int -> int -> int -> float -(** Get an integral using the Chemist's convention {% $(ij|kl)$ %}. *) - -val get_phys : t -> int -> int -> int -> int -> float -(** Get an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *) - val of_basis : Basis.t -> t (** Compute all ERI's for a given {!Basis.t}. *) -val to_file : ?cutoff:float -> filename:string -> t -> unit -(** Write the ERI's to a file, using the Physicist's convention. *) diff --git a/Nuclei/Nuclei.ml b/Nuclei/Nuclei.ml index 14acda5..1c4c14f 100644 --- a/Nuclei/Nuclei.ml +++ b/Nuclei/Nuclei.ml @@ -25,6 +25,12 @@ let of_xyz_file filename = +let of_zmt_string buffer = + Zmatrix.of_string buffer + |> Zmatrix.to_xyz + |> Array.map (fun (e,x,y,z) -> (e, Coordinate.angstrom_to_bohr (Angstrom.make {Point.x ; y ; z} ))) + + let of_zmt_file filename = let ic = open_in filename in let rec aux accu = @@ -36,9 +42,7 @@ let of_zmt_file filename = List.rev accu |> String.concat "\n" in aux [] - |> Zmatrix.of_string - |> Zmatrix.to_xyz - |> Array.map (fun (e,x,y,z) -> (e, Coordinate.angstrom_to_bohr (Angstrom.make {Point.x ; y ; z} ))) + |> of_zmt_string let to_string atoms = @@ -92,3 +96,30 @@ let charge nuclei = 0 nuclei |> Charge.of_int + +let to_xyz_string t = + [ string_of_int (Array.length t) ; "" ] @ + ( Array.mapi (fun i (e, coord) -> + let coord = + Coordinate.bohr_to_angstrom coord + in + Printf.sprintf " %5s %12.6f %12.6f %12.6f" + (Element.to_string e) coord.Angstrom.x coord.Angstrom.y coord.Angstrom.z + ) t + |> Array.to_list ) + |> String.concat "\n" + + +let to_t2_string t = + [ "# nAt nEl nCore nRyd" ; + Printf.sprintf " %d %d %d 0" (Array.length t) + (Array.fold_left (+) 0 (Array.map (fun (e,_) -> Element.to_int e) t) ) + (2 * Array.length t); + "# Znuc x y z" ] + @ (Array.mapi (fun i (e, coord) -> + Printf.sprintf " %5f %12.6f %12.6f %12.6f" + (Element.to_int e |> float_of_int) coord.Bohr.x coord.Bohr.y coord.Bohr.z + ) t + |> Array.to_list ) + |> String.concat "\n" + diff --git a/Nuclei/Nuclei.mli b/Nuclei/Nuclei.mli index 4aee0d5..abc9f59 100644 --- a/Nuclei/Nuclei.mli +++ b/Nuclei/Nuclei.mli @@ -8,6 +8,9 @@ type t = (Element.t * Bohr.t) array val of_xyz_file : string -> t (** Create from a file, in [xyz] format. *) +val of_zmt_string : string -> t +(** Create from a string, in z-matrix format. *) + val of_zmt_file : string -> t (** Create from a file, in z-matrix format. *) @@ -24,3 +27,5 @@ val repulsion : t -> float val charge : t -> Charge.t (** Sum of the charges of the nuclei. *) +val to_xyz_string : t -> string +val to_t2_string : t -> string diff --git a/SCF/Fock.ml b/SCF/Fock.ml index e801dca..ea69b93 100644 --- a/SCF/Fock.ml +++ b/SCF/Fock.ml @@ -14,6 +14,8 @@ type t = module Ao = AOBasis + + let make ~density ao_basis = let m_P = density and m_T = Lazy.force ao_basis.Ao.kin_ints |> KinInt.matrix @@ -27,6 +29,57 @@ let make ~density ao_basis = and m_J = Array.make_matrix nBas nBas 0. and m_K = Array.make_matrix nBas nBas 0. in +(* + let permutations i j k l = + let result = + [ [| i ; j ; k ; l |] ; + [| i ; l ; k ; j |] ; + [| k ; j ; i ; l |] ; + [| k ; l ; i ; j |] ; + [| j ; i ; l ; k |] ; + [| j ; k ; l ; i |] ; + [| l ; i ; j ; k |] ; + [| l ; k ; j ; i |] ] + in + if i<>k && j<>l && i<>j && k<>l then + result + else + List.map Zkey.of_int_array result + |> List.sort_uniq Pervasives.compare + |> List.map Zkey.to_int_array + in + let sum = ref 0 in + ERI.to_stream m_G + |> Stream.iter (fun { ERI.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } -> + permutations i_r1 j_r2 k_r1 l_r2 + |> List.iter ( fun ijkl -> + let mu = ijkl.(0) + and nu = ijkl.(2) + and lambda = ijkl.(1) + and sigma = ijkl.(3) + in + incr sum; + let p = m_P.{lambda,sigma} in + if abs_float p > epsilon then + let m_Jnu = m_J.(nu-1) in + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. p *. value ) + ); + + Printf.printf "%d %d\n%!" !sum (nBas*nBas*nBas*nBas); +*) + for sigma = 1 to nBas do + for nu = 1 to nBas do + let m_Jnu = m_J.(nu-1) in + for lambda = 1 to nBas do + let p = m_P.{lambda,sigma} in + for mu = 1 to nBas do + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. p *. + ERI.get_phys m_G mu lambda nu sigma + done + done + done + done; + (* for sigma = 1 to nBas do for nu = 1 to nBas do let m_Jnu = m_J.(nu-1) in @@ -45,6 +98,13 @@ let make ~density ao_basis = done done done; + for nu = 1 to nBas do + for mu = 1 to nu-1 do + m_J.(mu-1).(nu-1) <- m_J.(nu-1).(mu-1); + done + done; + *) + for nu = 1 to nBas do let m_Knu = m_K.(nu-1) in for sigma = 1 to nBas do @@ -58,12 +118,9 @@ let make ~density ao_basis = ERI.get_phys m_G mu lambda sigma nu done done - done - done; - for nu = 1 to nBas do + done; for mu = 1 to nu-1 do - m_J.(mu-1).(nu-1) <- m_J.(nu-1).(mu-1); - m_K.(mu-1).(nu-1) <- m_K.(nu-1).(mu-1); + m_K.(mu-1).(nu-1) <- m_Knu.(mu-1); done done; let m_J = Mat.of_array m_J diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index 2b7d6ad..de54bd8 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -173,7 +173,7 @@ let to_stream d = and k = ref 1 and l = ref 1 in - let f _ = + let f_dense _ = i := !i+1; if !i > !k then begin i := 1; @@ -195,7 +195,7 @@ let to_stream d = else None in - Stream.from f + Stream.from f_dense (** Write all integrals to a file with the convention *) diff --git a/Utils/FourIdxStorage.mli b/Utils/FourIdxStorage.mli index 065c31d..4736d02 100644 --- a/Utils/FourIdxStorage.mli +++ b/Utils/FourIdxStorage.mli @@ -25,12 +25,22 @@ val create : size:int -> [< `Dense | `Sparse ] -> t (** {2 Accessors} *) val get_chem : t -> int -> int -> int -> int -> float +(** Get an integral using the Chemist's convention {% $(ij|kl)$ %}. *) + val get_phys : t -> int -> int -> int -> int -> float +(** Get an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *) + val set_chem : t -> int -> int -> int -> int -> float -> unit +(** Set an integral using the Chemist's convention {% $(ij|kl)$ %}. *) + val set_phys : t -> int -> int -> int -> int -> float -> unit +(** Set an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *) val to_stream : t -> element Stream.t +(** Fetch all integrals from a stream. *) + (** {2 I/O} *) + val to_file : ?cutoff:float -> filename:string -> t -> unit (** Write the data to file, using the physicist's ordering. *) diff --git a/Utils/Zkey.ml b/Utils/Zkey.ml index 9e5f728..76a0820 100644 --- a/Utils/Zkey.ml +++ b/Utils/Zkey.ml @@ -31,6 +31,7 @@ let (<+) z x = type kind = | Three of Powers.t + | Four of (int * int * int * int) | Six of (Powers.t * Powers.t) | Nine of (Powers.t * Powers.t * Powers.t) | Twelve of (Powers.t * Powers.t * Powers.t * Powers.t) @@ -42,6 +43,13 @@ let of_powers_three { x=a ; y=b ; z=c ; _ } = ); make 3 a <+ b <+ c +let of_int_four i j k l = + assert ( + let alpha = i lor j lor k lor l in + alpha >= 0 && alpha < (1 lsl 15) + ); + make 4 i <+ j <+ k <+ l + let of_powers_six { x=a ; y=b ; z=c ; _ } { x=d ; y=e ; z=f ; _ } = assert ( let alpha = a lor b lor c lor d lor e lor f in @@ -76,12 +84,18 @@ let of_powers a = | Six (a,b) -> of_powers_six a b | Twelve (a,b,c,d) -> of_powers_twelve a b c d | Nine (a,b,c) -> of_powers_nine a b c + | _ -> invalid_arg "of_powers" let mask10 = 0x3ff and mask15 = 0x7fff +let of_int_array = function +| [| a ; b ; c ; d |] -> of_int_four a b c d +| _ -> invalid_arg "of_int_array" + + (** Transform the Zkey into an int array *) let to_int_array { left ; right ; kind } = match kind with @@ -91,6 +105,13 @@ let to_int_array { left ; right ; kind } = mask15 land right |] + | 4 -> [| + mask15 land (right lsr 45) ; + mask15 land (right lsr 30) ; + mask15 land (right lsr 15) ; + mask15 land right + |] + | 6 -> [| mask10 land (right lsr 50) ; mask10 land (right lsr 40) ; @@ -180,6 +201,7 @@ let to_powers { left ; right ; kind } = mask10 land (right lsr 10) , mask10 land right ) ) + | _ -> invalid_arg (__FILE__^": to_powers") diff --git a/Utils/Zkey.mli b/Utils/Zkey.mli index ec43a41..27475ad 100644 --- a/Utils/Zkey.mli +++ b/Utils/Zkey.mli @@ -33,6 +33,9 @@ val to_string : t -> string val of_powers_three : Powers.t -> t (** Create from a {!Powers.t}. *) +val of_int_four : int -> int -> int -> int -> t +(** Create from four integers. *) + val of_powers_six : Powers.t -> Powers.t -> t (** Create from two {!Powers.t}. *) @@ -44,6 +47,7 @@ val of_powers_twelve : Powers.t -> Powers.t -> Powers.t -> Powers.t -> t type kind = | Three of Powers.t + | Four of (int * int * int * int) | Six of (Powers.t * Powers.t) | Nine of (Powers.t * Powers.t * Powers.t) | Twelve of (Powers.t * Powers.t * Powers.t * Powers.t) @@ -54,6 +58,9 @@ val of_powers : kind -> t val to_int_array : t -> int array (** Convert to an int array. *) +val of_int_array : int array -> t +(** Convert from an int array. *) + val to_powers : t -> kind (** {1 Functions for hash tables} *)