From b72a1e251e77669af563e0b98ca3866b4aca4e97 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 4 Mar 2019 19:01:54 +0100 Subject: [PATCH] Added Hartree-Fock guess with projection --- Basis/AOBasis.ml | 4 ++-- Basis/KinInt.ml | 2 ++ Basis/MatrixOnBasis.mli | 3 +++ Basis/NucInt.ml | 3 +++ Basis/Overlap.ml | 44 +++++++++++++++++++++++++++++++++++++++++ MOBasis/MOBasis.ml | 41 +++++++++++++++++++++++++++++++++++++- MOBasis/MOBasis.mli | 6 ++++-- Makefile.include | 2 +- SCF/Guess.ml | 4 ++-- SCF/Guess.mli | 9 ++++++--- SCF/HartreeFock.ml | 26 ++++++++++++------------ SCF/HartreeFock.mli | 5 +++-- run_hartree_fock.ml | 42 ++++++++++++++++++++++++++++----------- 13 files changed, 154 insertions(+), 37 deletions(-) diff --git a/Basis/AOBasis.ml b/Basis/AOBasis.ml index 8f7ae13..5de5b24 100644 --- a/Basis/AOBasis.ml +++ b/Basis/AOBasis.ml @@ -7,8 +7,8 @@ type t = overlap : Overlap.t lazy_t; ortho : Orthonormalization.t lazy_t; eN_ints : NucInt.t lazy_t; - ee_ints : ERI.t lazy_t; kin_ints : KinInt.t lazy_t; + ee_ints : ERI.t lazy_t; cartesian : bool; } @@ -16,8 +16,8 @@ let basis t = t.basis let overlap t = Lazy.force t.overlap let ortho t = Lazy.force t.ortho let eN_ints t = Lazy.force t.eN_ints -let ee_ints t = Lazy.force t.ee_ints let kin_ints t = Lazy.force t.kin_ints +let ee_ints t = Lazy.force t.ee_ints let cartesian t = t.cartesian diff --git a/Basis/KinInt.ml b/Basis/KinInt.ml index bbb214e..8472721 100644 --- a/Basis/KinInt.ml +++ b/Basis/KinInt.ml @@ -161,6 +161,8 @@ let of_basis basis = Mat.detri result; result +let of_basis_pair first_basis second_basis = + failwith "Not implemented" (** Write all kinetic integrals to a file *) let to_file ~filename kinetic = diff --git a/Basis/MatrixOnBasis.mli b/Basis/MatrixOnBasis.mli index 642d362..99a8f53 100644 --- a/Basis/MatrixOnBasis.mli +++ b/Basis/MatrixOnBasis.mli @@ -7,6 +7,9 @@ type t val of_basis : Basis.t -> t (** Build from a {!Basis.t}. *) +val of_basis_pair : Basis.t -> Basis.t -> t +(** Build from a pair of {!Basis.t}. *) + val to_file : filename:string -> t -> unit (** Write the integrals in a file. *) diff --git a/Basis/NucInt.ml b/Basis/NucInt.ml index 86d8a9a..2001d18 100644 --- a/Basis/NucInt.ml +++ b/Basis/NucInt.ml @@ -114,3 +114,6 @@ let to_file ~filename eni_array = let of_basis basis = invalid_arg "of_basis_nuclei should be called for NucInt" +let of_basis_pair basis = + failwith "Not implemented" + diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index 756fa17..f73dfa0 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -130,6 +130,50 @@ let of_basis basis = result +(** Create mixed overlap matrix *) +let of_basis_pair first_basis second_basis = + let to_powers x = + let open Zkey in + match to_powers x with + | Three x -> x + | _ -> assert false + in + + let n = Bs.size first_basis + and m = Bs.size second_basis + and first = Bs.contracted_shells first_basis + and second = Bs.contracted_shells second_basis + in + + let result = Mat.create n m in + for j=0 to (Array.length second) - 1 do + for i=0 to (Array.length first) - 1 do + (* Compute all the integrals of the class *) + let cls = + contracted_class first.(i) second.(j) + in + + Array.iteri (fun j_c powers_j -> + let j_c = Cs.index second.(j) + j_c + 1 in + let xj = to_powers powers_j in + Array.iteri (fun i_c powers_i -> + let i_c = Cs.index first.(i) + i_c + 1 in + let xi = to_powers powers_i in + let key = + Zkey.of_powers_six xi xj + in + let value = + try Zmap.find cls key + with Not_found -> 0. + in + result.{i_c,j_c} <- value; + ) (Am.zkey_array (Singlet (Cs.ang_mom first.(i)))) + ) (Am.zkey_array (Singlet (Cs.ang_mom second.(j)))) + done; + done; + result + + (** Write all overlap integrals to a file *) let to_file ~filename overlap = diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml index 40ec33c..d32d457 100644 --- a/MOBasis/MOBasis.ml +++ b/MOBasis/MOBasis.ml @@ -8,7 +8,7 @@ module HF = HartreeFock module Si = Simulation type mo_type = - | RHF | ROHF | UHF | CASSCF + | RHF | ROHF | UHF | CASSCF | Projected | Natural of string | Localized of string @@ -200,6 +200,45 @@ let of_hartree_fock hf = make ~simulation ~mo_type ~mo_occupation ~mo_coef () +let of_mo_basis simulation other = + + let mo_coef = + let basis = Simulation.ao_basis simulation in + let basis_other = ao_basis other in + let m_S = + Overlap.(matrix @@ of_basis_pair + (AOBasis.basis basis) + (AOBasis.basis basis_other) ) + in + let m_X = AOBasis.ortho basis in + (* Project other vectors in the current basis *) + let m_C = + gemm m_S @@ mo_coef other + in + (* Append dummy vectors to the input vectors *) + let result = + let vecs = Mat.to_col_vecs m_X in + Array.iteri (fun i v -> if (i < Array.length vecs) then vecs.(i) <- v) + (Mat.to_col_vecs m_C) ; + Mat.of_col_vecs vecs + in + (* Gram-Schmidt Orthonormalization *) + gemm m_X @@ (Util.qr_ortho @@ gemm ~transa:`T m_X result) + in + + let mo_occupation = + let occ = mo_occupation other in + Vec.init (Mat.dim2 mo_coef) (fun i -> + if (i <= Vec.dim occ) then occ.{i} + else 0.) + in + make ~simulation ~mo_type:Projected ~mo_occupation ~mo_coef () + + + + + + let pp_mo ?(start=1) ?finish ppf t = let open Lacaml.Io in diff --git a/MOBasis/MOBasis.mli b/MOBasis/MOBasis.mli index 49d7890..124e531 100644 --- a/MOBasis/MOBasis.mli +++ b/MOBasis/MOBasis.mli @@ -7,7 +7,7 @@ open Lacaml.D type mo_type = - | RHF | ROHF | UHF | CASSCF + | RHF | ROHF | UHF | CASSCF | Projected | Natural of string | Localized of string @@ -51,7 +51,7 @@ val size : t -> int (** {1 Creators} *) -val make : simulation:Simulation.t -> +val make : simulation:Simulation.t -> mo_type:mo_type -> mo_occupation:Vec.t -> mo_coef:Mat.t -> @@ -61,6 +61,8 @@ val make : simulation:Simulation.t -> val of_hartree_fock : HartreeFock.t -> t (** Build MOs from a Restricted Hartree-Fock calculation. *) +val of_mo_basis : Simulation.t -> t -> t +(** Project the MOs of the other basis on the current one. *) (** {1 Printers} *) diff --git a/Makefile.include b/Makefile.include index ab5345d..4adbe57 100644 --- a/Makefile.include +++ b/Makefile.include @@ -1,6 +1,6 @@ .NOPARALLEL: -INCLUDE_DIRS=Parallel,Nuclei,Utils,Basis,SCF,MOBasis,CI +INCLUDE_DIRS=Parallel,Nuclei,Utils,Basis,SCF,MOBasis,CI,F12 LIBS= PKGS= OCAMLBUILD=ocamlbuild -j 0 -cflags $(ocamlcflags) -lflags $(ocamllflags) $(ocamldocflags) -Is $(INCLUDE_DIRS) -ocamlopt $(ocamloptflags) $(mpi) diff --git a/SCF/Guess.ml b/SCF/Guess.ml index e6d81fa..aa8d1e4 100644 --- a/SCF/Guess.ml +++ b/SCF/Guess.ml @@ -4,6 +4,7 @@ open Util type guess = | Hcore of Mat.t | Huckel of Mat.t +| Matrix of Mat.t type t = guess @@ -49,8 +50,7 @@ let make ?(nocc=0) ~guess ao_basis = match guess with | `Hcore -> Hcore (hcore_guess ao_basis) | `Huckel -> Huckel (huckel_guess ao_basis nocc) - - + | `Matrix m -> Matrix m diff --git a/SCF/Guess.mli b/SCF/Guess.mli index 25fdfed..3961adc 100644 --- a/SCF/Guess.mli +++ b/SCF/Guess.mli @@ -1,13 +1,16 @@ +open Lacaml.D + (** Guess for Hartree-Fock calculations. *) type guess = -| Hcore of Lacaml.D.Mat.t -| Huckel of Lacaml.D.Mat.t +| Hcore of Mat.t (* Core Hamiltonian Matrix *) +| Huckel of Mat.t (* Huckel Hamiltonian Matrix *) +| Matrix of Mat.t (* Guess Eigenvectors *) type t = guess -val make : ?nocc:int -> guess:[ `Hcore | `Huckel ] -> AOBasis.t -> t +val make : ?nocc:int -> guess:[ `Hcore | `Huckel | `Matrix of Mat.t ] -> AOBasis.t -> t (** {2 Tests} *) diff --git a/SCF/HartreeFock.ml b/SCF/HartreeFock.ml index 809b5d6..273466b 100644 --- a/SCF/HartreeFock.ml +++ b/SCF/HartreeFock.ml @@ -211,10 +211,6 @@ let make Si.ao_basis simulation in - (* Initial guess *) - let guess = - Guess.make ~nocc ~guess ao_basis - in (* Orthogonalization matrix *) let m_X = @@ -236,17 +232,21 @@ let make in (* Guess coefficients *) - let m_C = - let m_H = - match guess with - | Guess.Hcore m_H -> m_H - | Guess.Huckel m_H -> m_H - in - let m_Hmo = xt_o_x m_H m_X in - let m_C', _ = diagonalize_symm m_Hmo in - gemm m_X m_C' + let guess = + Guess.make ~nocc ~guess ao_basis in + let m_C = + let c_of_h m_H = + let m_Hmo = xt_o_x m_H m_X in + let m_C', _ = diagonalize_symm m_Hmo in + gemm m_X m_C' + in + match guess with + | Guess.Hcore m_H -> c_of_h m_H + | Guess.Huckel m_H -> c_of_h m_H + | Guess.Matrix m_C -> m_C + in (* A single SCF iteration *) let scf_iteration_rhf data = diff --git a/SCF/HartreeFock.mli b/SCF/HartreeFock.mli index 8012eeb..3384e2a 100644 --- a/SCF/HartreeFock.mli +++ b/SCF/HartreeFock.mli @@ -57,9 +57,10 @@ val empty: hartree_fock_data val make : ?kind:hartree_fock_kind -> - ?guess:[ `Hcore | `Huckel ] -> + ?guess:[ `Hcore | `Huckel | `Matrix of Mat.t ] -> ?max_scf:int -> - ?level_shift:float -> ?threshold_SCF:float -> Simulation.t -> t + ?level_shift:float -> ?threshold_SCF:float -> + Simulation.t -> t (** {1 Printers} *) diff --git a/run_hartree_fock.ml b/run_hartree_fock.ml index 4a81379..c985076 100644 --- a/run_hartree_fock.ml +++ b/run_hartree_fock.ml @@ -1,3 +1,5 @@ +open Lacaml.D + let () = let open Command_line in begin @@ -19,10 +21,18 @@ let () = { short='c' ; long="charge" ; opt=Optional; arg=With_arg ""; doc="Total charge of the molecule. Default is 0"; } ; + + { short='g' ; long="guess" ; opt=Optional; + arg=With_arg ""; + doc="Guess with another basis set."; } ; ] end; - (* Handle options *) + let ppf = + if Parallel.master then Format.std_formatter + else Printing.ppf_dev_null + in + let basis_file = Util.of_some @@ Command_line.get "basis" in let nuclei_file = Util.of_some @@ Command_line.get "xyz" in @@ -39,20 +49,30 @@ let () = | None -> 1 in - let s = - Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file + let s = Simulation.of_filenames + ~charge ~multiplicity ~nuclei:nuclei_file basis_file in + + let guess = + match Command_line.get "guess" with + | None -> `Huckel + | Some filename -> + let s_guess = Simulation.of_filenames + ~charge ~multiplicity ~nuclei:nuclei_file filename in + let hf = HartreeFock.make s_guess in + Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf; + let guess_mos = + MOBasis.of_hartree_fock hf + |> MOBasis.of_mo_basis s + in + `Matrix (MOBasis.mo_coef guess_mos) in - let hf = - HartreeFock.make s - in + let hf = HartreeFock.make ~guess s in - let ppf = - if Parallel.master then Format.std_formatter - else Printing.ppf_dev_null - in - Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf; + Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf + (* let mos = MOBasis.of_hartree_fock hf in Format.fprintf ppf "@[%a@]@." (fun ppf x -> MOBasis.pp_mo ppf x) mos + *)