From 4f41af9e31e318b1f060586f0ae9d61026dbdcaa Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 20 Jul 2018 16:09:06 +0200 Subject: [PATCH] Four-idx transformation --- Basis/KinInt.ml | 1 + Basis/MatrixOnBasis.mli | 5 +- Basis/NucInt.ml | 1 + Basis/Overlap.ml | 1 + MOBasis/MOBasis.ml | 181 +++++++++++++++++++++++++++++++++++++++ MOBasis/MOBasis.mli | 43 ++++++++++ Makefile.include | 2 +- SCF/HartreeFock_type.ml | 33 ++++--- SCF/HartreeFock_type.mli | 9 +- SCF/RHF.ml | 6 +- Utils/FourIdxStorage.ml | 17 ++++ Utils/FourIdxStorage.mli | 12 +++ Utils/Util.ml | 6 +- Utils/Util.mli | 7 +- 14 files changed, 307 insertions(+), 17 deletions(-) create mode 100644 MOBasis/MOBasis.ml create mode 100644 MOBasis/MOBasis.mli diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index daff5f8..bbb214e 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -13,6 +13,7 @@ module Psp = PrimitiveShellPair type t = Mat.t external matrix : t -> Mat.t = "%identity" +external of_matrix : Mat.t -> t = "%identity" let cutoff = integrals_cutoff diff --git a/Basis/MatrixOnBasis.mli b/Basis/MatrixOnBasis.mli index c8955e8..642d362 100644 --- a/Basis/MatrixOnBasis.mli +++ b/Basis/MatrixOnBasis.mli @@ -11,4 +11,7 @@ val to_file : filename:string -> t -> unit (** Write the integrals in a file. *) val matrix : t -> Mat.t -(** Returns the matrix suitable for Lacaml *) +(** Returns the matrix suitable for Lacaml. *) + +val of_matrix : Mat.t -> t +(** Build from a Lacaml matrix. *) diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml index c5baff4..86d8a9a 100644 --- a/Basis/NucInt.ml +++ b/Basis/NucInt.ml @@ -6,6 +6,7 @@ open Lacaml.D type t = Mat.t external matrix : t -> Mat.t = "%identity" +external of_matrix : Mat.t -> t = "%identity" module Am = AngularMomentum module Bs = Basis diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index 7451cbf..756fa17 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -4,6 +4,7 @@ open Lacaml.D type t = Mat.t external matrix : t -> Mat.t = "%identity" +external of_matrix : Mat.t -> t = "%identity" module Am = AngularMomentum diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml new file mode 100644 index 0000000..c3c3db9 --- /dev/null +++ b/MOBasis/MOBasis.ml @@ -0,0 +1,181 @@ +open Lacaml.D +open Util + +(** One-electron orthogonal basis set, corresponding to Molecular Orbitals. *) + +module HF = HartreeFock_type +module Si = Simulation + +type mo_class = +| Core of int +| Inactive of int +| Active of int +| Virtual of int +| Deleted of int + +type mo_type = +| RHF | ROHF | CASSCF +| Natural of string +| Localized of string + +type t = +{ + ao_basis : AOBasis.t; (* Atomic basis set on which the MOs are built. *) + mo_type : mo_type; (* Kind of MOs (RHF, CASSCF, Localized... *) + mo_class : mo_class array; (* CI-Class of the MOs *) + mo_occupation : Vec.t; (* Occupation numbers *) + mo_coef : Mat.t; (* Matrix of the MO coefficients in the AO basis *) + eN_ints : NucInt.t lazy_t; (* Electron-nucleus potential integrals *) + ee_ints : ERI.t lazy_t; (* Electron-electron potential integrals *) + kin_ints : KinInt.t lazy_t; (* Kinetic energy integrals *) +} + + + +let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = + gemm ~transa:`T mo_coef @@ + gemm ao_matrix mo_coef + + +let ao_matrix_of_mo_matrix ~mo_coef ~ao_overlap mo_matrix = + let sc = gemm ao_overlap mo_coef in + gemm sc @@ + gemm ~transb:`T mo_matrix sc + + +let four_index_transform ~mo_coef eri_ao = + + let ao_num = Mat.dim1 mo_coef in + let mo_num = Mat.dim2 mo_coef in + let eri_mo = ERI.create ~size:mo_num `Dense in + + let mo_num_2 = mo_num * mo_num in + let ao_num_2 = ao_num * ao_num in + let ao_mo_num = ao_num * mo_num in + + let range_mo = list_range ~start:1 mo_num in + let range_ao = list_range ~start:1 ao_num in + + let u = + Mat.create mo_num_2 mo_num + and o = + Mat.create ao_num ao_num_2 + and p = + Mat.create ao_num_2 mo_num + and q = + Mat.create ao_mo_num mo_num + in + Printf.eprintf "Transforming %d integrals : %!" mo_num; + List.iter (fun delta -> + Printf.eprintf "%d %!" delta; + Mat.fill u 0.; + + List.iter (fun l -> + + let jk = ref 0 in + List.iter (fun k -> + List.iter (fun j -> + incr jk; + ERI.get_phys_all_i eri_ao ~j ~k ~l + |> Array.iteri (fun i x -> o.{i+1,!jk} <- x) + ) range_ao + ) range_ao; + (* o_i_jk *) + + let p = + gemm ~transa:`T ~c:p o mo_coef + (* p_jk_alpha = \sum_i o_i_jk c_i_alpha *) + in + let p' = + Bigarray.reshape_2 (Bigarray.genarray_of_array2 p) ao_num ao_mo_num + (* p_j_kalpha *) + in + + let q = + gemm ~transa:`T ~c:q p' mo_coef + (* q_kalpha_beta = \sum_j p_j_kalpha c_j_beta *) + in + let q' = + Bigarray.reshape_2 (Bigarray.genarray_of_array2 q) ao_num mo_num_2 + (* q_k_alphabeta = \sum_j p_j_kalpha c_j_beta *) + in + + ignore @@ + gemm ~transa:`T ~beta:1. ~c:u q' mo_coef + (* u_alphabeta_gamma = \sum_k q_k_alphabeta c_k_gamma *) + + ) range_ao; + let u = + Bigarray.reshape + (Bigarray.genarray_of_array2 u) + [| mo_num ; mo_num ; mo_num |] + |> Bigarray.array3_of_genarray + in + List.iter (fun gamma -> + List.iter (fun beta -> + List.iter (fun alpha -> + let x = u.{alpha,beta,gamma} in + if x <> 0. then + ERI.set_phys eri_mo alpha beta gamma delta x + ) range_mo + ) range_mo + ) range_mo + ) range_mo; + Printf.eprintf "\n%!"; + + eri_mo + + +let make ~ao_basis ~mo_type ~mo_class ~mo_occupation ~mo_coef () = + let eN_ints = lazy ( + Lazy.force ao_basis.AOBasis.eN_ints + |> NucInt.matrix + |> mo_matrix_of_ao_matrix ~mo_coef + |> NucInt.of_matrix + ) + and kin_ints = lazy ( + Lazy.force ao_basis.AOBasis.kin_ints + |> KinInt.matrix + |> mo_matrix_of_ao_matrix ~mo_coef + |> KinInt.of_matrix + ) + and ee_ints = lazy ( + Lazy.force ao_basis.AOBasis.ee_ints + |> four_index_transform ~mo_coef + ) + in + { ao_basis ; mo_type ; mo_class ; mo_occupation ; mo_coef ; + eN_ints ; ee_ints ; kin_ints } + + +let of_rhf ~frozen_core hf = + let simulation = hf.HF.simulation in + let nocc = hf.HF.nocc in + let ncore = + if frozen_core then Nuclei.small_core simulation.Si.nuclei + else 0 + in + let mo_num = Vec.dim hf.HF.eigenvalues in + + let ao_basis = simulation.Si.ao_basis in + let mo_type = RHF in + let mo_class = + Array.init mo_num (fun i -> + if (i < ncore) then Core i + else + if (i < nocc ) then Inactive i + else Virtual i) + in + let mo_occupation = + Array.init mo_num (fun i -> + if i < nocc then 2. else 0.) + |> Vec.of_array + in + let mo_coef = hf.HF.eigenvectors in + make ~ao_basis ~mo_type ~mo_class ~mo_occupation ~mo_coef () + + +let of_hartree_fock ~frozen_core = function +| HF.RHF hf -> of_rhf ~frozen_core hf +| _ -> assert false + diff --git a/MOBasis/MOBasis.mli b/MOBasis/MOBasis.mli new file mode 100644 index 0000000..07e5f76 --- /dev/null +++ b/MOBasis/MOBasis.mli @@ -0,0 +1,43 @@ +(** Data structure to represent the molecular orbitals. *) +open Lacaml.D + +type mo_class = +| Core of int +| Inactive of int +| Active of int +| Virtual of int +| Deleted of int + + +type mo_type = +| RHF | ROHF | CASSCF +| Natural of string +| Localized of string + + +type t = private +{ + ao_basis : AOBasis.t; (* Atomic basis set on which the MOs are built. *) + mo_type : mo_type; (* Kind of MOs (RHF, CASSCF, Localized... *) + mo_class : mo_class array; (* CI-Class of the MOs *) + mo_occupation : Vec.t; (* Occupation numbers *) + mo_coef : Mat.t; (* Matrix of the MO coefficients in the AO basis *) + eN_ints : NucInt.t lazy_t; (* Electron-nucleus potential integrals *) + ee_ints : ERI.t lazy_t; (* Electron-electron potential integrals *) + kin_ints : KinInt.t lazy_t; (* Kinetic energy integrals *) +} + + + +val make : ao_basis:AOBasis.t -> + mo_type:mo_type -> + mo_class:mo_class array -> + mo_occupation:Vec.t -> + mo_coef:Mat.t -> + unit -> t +(** Function to build a data structure representing the molecular orbitals. *) + +val of_hartree_fock : frozen_core:bool -> HartreeFock_type.t -> t +(** Build MOs from a Restricted Hartree-Fock calculation. *) + + diff --git a/Makefile.include b/Makefile.include index afec22d..2c79cc4 100644 --- a/Makefile.include +++ b/Makefile.include @@ -1,6 +1,6 @@ .NOPARALLEL: -INCLUDE_DIRS=Nuclei,Utils,Basis,SCF +INCLUDE_DIRS=Nuclei,Utils,Basis,SCF,MOBasis LIBS= PKGS= OCAMLBUILD=ocamlbuild -j 0 -cflags $(ocamlcflags) -lflags $(ocamlcflags) $(ocamldocflags) -Is $(INCLUDE_DIRS) -ocamlopt $(ocamloptflags) diff --git a/SCF/HartreeFock_type.ml b/SCF/HartreeFock_type.ml index 0943ef0..ca19393 100644 --- a/SCF/HartreeFock_type.ml +++ b/SCF/HartreeFock_type.ml @@ -1,8 +1,9 @@ open Lacaml.D open Util -type t = +type s = { + simulation : Simulation.t; guess : Guess.t; eigenvectors : Lacaml.D.Mat.t ; eigenvalues : Lacaml.D.Vec.t ; @@ -17,6 +18,12 @@ type t = } +type t = +| RHF of s (** Restricted Hartree-Fock *) +| ROHF of s (** Restricted Open-shell Hartree-Fock *) +| UHF of s (** Unrestricted Hartree-Fock *) + + let iterations_to_string hf_calc = " # HF energy Convergence HOMO-LUMO ---------------------------------------------------" :: @@ -100,15 +107,21 @@ let mos_to_string hf_calc = " ] -let to_string hf_calc = - String.concat "\n" [ " - ===================================================== - Hartree-Fock - =====================================================" ; "" ; - iterations_to_string hf_calc ; "" ; - summary hf_calc ; "" ; - mos_to_string hf_calc ; "" ; - ] +let to_string hf = + let aux hf_calc r = + String.concat "\n" [ Printf.sprintf " + ===================================================== + %s Hartree-Fock + =====================================================" r ; "" ; + iterations_to_string hf_calc ; "" ; + summary hf_calc ; "" ; + mos_to_string hf_calc ; "" ; + ] + in + match hf with + | RHF hf_calc -> aux hf_calc "Restricted" + | UHF hf_calc -> aux hf_calc "Unrestricted" + | ROHF hf_calc -> aux hf_calc "Restricted Open-shell" diff --git a/SCF/HartreeFock_type.mli b/SCF/HartreeFock_type.mli index 3d4fa01..a79d046 100644 --- a/SCF/HartreeFock_type.mli +++ b/SCF/HartreeFock_type.mli @@ -1,7 +1,8 @@ (** Data structure representing the output of a Hartree-Fock caculation *) -type t = +type s = { + simulation : Simulation.t; (** Simulation which was used for HF calculation *) guess : Guess.t; (** Initial guess *) eigenvectors : Lacaml.D.Mat.t ; (** Final eigenvectors *) eigenvalues : Lacaml.D.Vec.t ; (** Final eigenvalues *) @@ -16,6 +17,12 @@ type t = (** Energy, convergence and HOMO-LUMO gap of all iterations *) } +type t = +| RHF of s (** Restricted Hartree-Fock *) +| ROHF of s (** Restricted Open-shell Hartree-Fock *) +| UHF of s (** Unrestricted Hartree-Fock *) + + val to_string : t -> string (** Results of a Hartree-Fock calculation pretty-printed in a string. *) diff --git a/SCF/RHF.ml b/SCF/RHF.ml index 0deeca4..f29f0f7 100644 --- a/SCF/RHF.ml +++ b/SCF/RHF.ml @@ -158,7 +158,9 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift= List.rev ( (energy, error, gap) :: iterations ) |> Array.of_list in - { HartreeFock_type. + HartreeFock_type.(RHF + { + simulation; nocc; guess ; eigenvectors = m_C ; @@ -170,7 +172,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift= eN_energy = Mat.gemm_trace m_P m_V; coulomb_energy = 0.5 *. Mat.gemm_trace m_P m_J; exchange_energy = 0.5 *. Mat.gemm_trace m_P m_K; - } + }) in diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index c5e8da6..76532cc 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -166,6 +166,23 @@ type element = (** Element for the stream *) } +let get_phys_all_i d ~j ~k ~l = + Array.init d.size (fun i -> get_phys d (i+1) j k l) + + +let get_chem_all_i d ~j ~k ~l = + Array.init d.size (fun i -> get_chem d (i+1) j k l) + + +let get_phys_all_ji d ~k ~l = + Array.init d.size (fun j -> get_phys_all_i d ~j:(j+1) ~k ~l) + + +let get_chem_all_ji d ~k ~l = + Array.init d.size (fun j -> get_chem_all_i d ~j:(j+1) ~k ~l) + + + let to_stream d = let i = ref 0 diff --git a/Utils/FourIdxStorage.mli b/Utils/FourIdxStorage.mli index 35eef2f..3129200 100644 --- a/Utils/FourIdxStorage.mli +++ b/Utils/FourIdxStorage.mli @@ -36,6 +36,18 @@ val set_chem : t -> int -> int -> int -> int -> float -> unit val set_phys : t -> int -> int -> int -> int -> float -> unit (** Set an integral using the Physicist's convention {% $\langle ij|kl \rangle$ %}. *) +val get_chem_all_i : t -> j:int -> k:int -> l:int -> float array +(** Get all integrals in an array [a.(i-1) =] {% $(\cdot j|kl)$ %} . *) + +val get_phys_all_i : t -> j:int -> k:int -> l:int -> float array +(** Get all integrals in an array [a.(i-1) =] {% $\langle \cdot j|kl \rangle$ %} . *) + +val get_chem_all_ji : t -> k:int -> l:int -> float array array +(** Get all integrals in an array [a.(j-1).(i-1) =] {% $(\cdot \cdot|kl)$ %} . *) + +val get_phys_all_ji : t -> k:int -> l:int -> float array array +(** Get all integrals in an array [a.(j-1).(i-1) =] {% $\langle \cdot \cdot|kl \rangle$ %} . *) + val to_stream : t -> element Stream.t (** Retrun the data structure as a stream. *) diff --git a/Utils/Util.ml b/Utils/Util.ml index ca52ead..97168fe 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -169,7 +169,7 @@ let list_some l = (** {2 Stream functions} *) -let range ?(start=0) n = +let stream_range ?(start=0) n = Stream.from (fun i -> let result = i+start in if result < n then @@ -178,6 +178,10 @@ let range ?(start=0) n = ) +let list_range ?(start=0) n = + Array.init n (fun i -> start+i) |> Array.to_list + + (** {2 Linear algebra} *) diff --git a/Utils/Util.mli b/Utils/Util.mli index 663638a..4314561 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -63,8 +63,13 @@ val list_some : 'a option list -> 'a list (** Filters out all [None] elements of the list, and returns the elements without the [Some]. *) +val list_range : ?start:int -> int -> int list +(** Returns a list [start ; start+1 ; ... ; start+(n-1)]. Default is [start=0]. *) + (** {2 Useful streams} *) -val range : ?start:int -> int -> int Stream.t +val stream_range : ?start:int -> int -> int Stream.t +(** Returns a stream . Default is [start=0]. *) + (** {2 Linear algebra } *)