From 3482e02695eb4d4db37bae4fd5ad5818c66aa8c0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 31 May 2018 16:46:45 +0200 Subject: [PATCH] Working on DIIS --- .merlin | 2 +- HartreeFock/Guess.ml | 17 ------ HartreeFock/Guess.mli | 12 ----- Makefile.include | 2 +- {HartreeFock => SCF}/Fock.ml | 0 SCF/Guess.ml | 49 +++++++++++++++++ SCF/Guess.mli | 13 +++++ {HartreeFock => SCF}/HartreeFock.ml | 0 {HartreeFock => SCF}/HartreeFock_type.ml | 0 {HartreeFock => SCF}/HartreeFock_type.mli | 0 {HartreeFock => SCF}/RHF.ml | 64 ++++++++++++++++------- Simulation.ml | 2 +- Utils/DIIS.ml | 60 +++++++++++++++++++++ Utils/DIIS.mli | 16 +++--- Utils/Util.ml | 29 ++++++++-- 15 files changed, 205 insertions(+), 61 deletions(-) delete mode 100644 HartreeFock/Guess.ml delete mode 100644 HartreeFock/Guess.mli rename {HartreeFock => SCF}/Fock.ml (100%) create mode 100644 SCF/Guess.ml create mode 100644 SCF/Guess.mli rename {HartreeFock => SCF}/HartreeFock.ml (100%) rename {HartreeFock => SCF}/HartreeFock_type.ml (100%) rename {HartreeFock => SCF}/HartreeFock_type.mli (100%) rename {HartreeFock => SCF}/RHF.ml (69%) create mode 100644 Utils/DIIS.ml diff --git a/.merlin b/.merlin index 05d34eb..f2649f0 100644 --- a/.merlin +++ b/.merlin @@ -3,5 +3,5 @@ S . S Nuclei S Utils S Basis -S HartreeFock +S SCF B _build/** diff --git a/HartreeFock/Guess.ml b/HartreeFock/Guess.ml deleted file mode 100644 index 8a4c471..0000000 --- a/HartreeFock/Guess.ml +++ /dev/null @@ -1,17 +0,0 @@ -open Lacaml.D - -type guess = -| Hcore of Mat.t - -type t = guess - - -let make ?guess:(guess=`Hcore) simulation = - let eN_ints = Lazy.force simulation.Simulation.eN_ints - and kin_ints = Lazy.force simulation.Simulation.kin_ints - in - match guess with - | `Hcore -> Hcore (Mat.add eN_ints kin_ints) - - - diff --git a/HartreeFock/Guess.mli b/HartreeFock/Guess.mli deleted file mode 100644 index 499822d..0000000 --- a/HartreeFock/Guess.mli +++ /dev/null @@ -1,12 +0,0 @@ -(** Guess for Hartree-Fock calculations. *) - -type guess = -| Hcore of Lacaml.D.Mat.t - -type t = guess - - -val make : ?guess:[ `Hcore ] -> Simulation.t -> t - - - diff --git a/Makefile.include b/Makefile.include index 99a6ed9..cfce256 100644 --- a/Makefile.include +++ b/Makefile.include @@ -1,6 +1,6 @@ .NOPARALLEL: -INCLUDE_DIRS=Nuclei,Utils,Basis,HartreeFock +INCLUDE_DIRS=Nuclei,Utils,Basis,SCF LIBS= PKGS= OCAMLBUILD=ocamlbuild -j 0 -cflags $(ocamlcflags) -lflags $(ocamlcflags) $(ocamldocflags) -Is $(INCLUDE_DIRS) -ocamlopt $(ocamloptflags) diff --git a/HartreeFock/Fock.ml b/SCF/Fock.ml similarity index 100% rename from HartreeFock/Fock.ml rename to SCF/Fock.ml diff --git a/SCF/Guess.ml b/SCF/Guess.ml new file mode 100644 index 0000000..653f193 --- /dev/null +++ b/SCF/Guess.ml @@ -0,0 +1,49 @@ +open Lacaml.D +open Util + +type guess = +| Hcore of Mat.t +| Huckel of Mat.t + +type t = guess + +module Si = Simulation +module El = Electrons + +let hcore_guess simulation = + let eN_ints = Lazy.force simulation.Si.eN_ints + and kin_ints = Lazy.force simulation.Si.kin_ints + in + Mat.add eN_ints kin_ints + + +let huckel_guess simulation = + let c = 0.5 *. 1.75 in + let ao_num = Basis.size simulation.Si.basis in + let eN_ints = Lazy.force simulation.Si.eN_ints + and kin_ints = Lazy.force simulation.Si.kin_ints + and overlap = Lazy.force simulation.Si.overlap + and m_X = Lazy.force simulation.Si.overlap_ortho + in + let diag = Array.init (ao_num+1) (fun i -> if i=0 then 0. else + eN_ints.{i,i} +. kin_ints.{i,i}) + in + + let nocc = + simulation.Si.electrons.El.n_alpha + in + let m_F = (Fock.make ~density:(gemm ~alpha:2. ~transb:`T ~k:nocc m_X m_X) simulation).Fock.fock in + for j=1 to ao_num do + for i=1 to ao_num do + if (i <> j) then + m_F.{i,j} <- c *. overlap.{i,j} *. (diag.(i) +. diag.(j)) (*TODO Pseudo *) + done; + done; + m_F + +let make ~guess simulation = + match guess with + | `Hcore -> Hcore (hcore_guess simulation) + | `Huckel -> Huckel (huckel_guess simulation) + + diff --git a/SCF/Guess.mli b/SCF/Guess.mli new file mode 100644 index 0000000..2fa54c2 --- /dev/null +++ b/SCF/Guess.mli @@ -0,0 +1,13 @@ +(** Guess for Hartree-Fock calculations. *) + +type guess = +| Hcore of Lacaml.D.Mat.t +| Huckel of Lacaml.D.Mat.t + +type t = guess + + +val make : guess:[ `Hcore | `Huckel ] -> Simulation.t -> t + + + diff --git a/HartreeFock/HartreeFock.ml b/SCF/HartreeFock.ml similarity index 100% rename from HartreeFock/HartreeFock.ml rename to SCF/HartreeFock.ml diff --git a/HartreeFock/HartreeFock_type.ml b/SCF/HartreeFock_type.ml similarity index 100% rename from HartreeFock/HartreeFock_type.ml rename to SCF/HartreeFock_type.ml diff --git a/HartreeFock/HartreeFock_type.mli b/SCF/HartreeFock_type.mli similarity index 100% rename from HartreeFock/HartreeFock_type.mli rename to SCF/HartreeFock_type.mli diff --git a/HartreeFock/RHF.ml b/SCF/RHF.ml similarity index 69% rename from HartreeFock/RHF.ml rename to SCF/RHF.ml index e072328..c6a778d 100644 --- a/HartreeFock/RHF.ml +++ b/SCF/RHF.ml @@ -2,7 +2,7 @@ open Util open Lacaml.D open Simulation -let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64) +let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=0.1) ?threshold_SCF:(threshold_SCF=1.e-6) simulation = (* Number of occupied MOs *) @@ -21,7 +21,7 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64) (* Orthogonalization matrix *) let m_X = - Lazy.force simulation.overlap_ortho + Lazy.force simulation.overlap_ortho in @@ -34,8 +34,16 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64) and m_V = Lazy.force simulation.eN_ints in + (* Level shift *) + let m_LS = + Array.init (Mat.dim2 m_X) (fun i -> + if i > nocc then level_shift else 0.) + |> Vec.of_array + |> Mat.of_diag + in + (* SCF iterations *) - let rec loop nSCF iterations m_C = + let rec loop nSCF iterations m_C diis = (* Density matrix over nocc occupied MOs *) let m_P = @@ -50,9 +58,31 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64) x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange in + let error_fock = + let fps = + gemm m_F (gemm m_P m_S) + and spf = + gemm m_S (gemm m_P m_F) + in + Mat.sub fps spf + in + + let diis = + DIIS.append ~p:(Mat.as_vec m_F) ~e:(Mat.as_vec error_fock) diis + in + + let m_F_diis = + let x = + Bigarray.genarray_of_array1 (DIIS.next diis) + in + Bigarray.reshape_2 x (Mat.dim1 m_F) (Mat.dim2 m_F) + in + + (* Fock matrix in MO basis *) let m_Fmo = - xt_o_x m_F m_X + xt_o_x m_F_diis m_C + |> Mat.add m_LS in (* MOs in old MO basis *) @@ -62,7 +92,7 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64) (* MOs in AO basis *) let m_C = - gemm m_X m_C' + gemm m_C m_C' in (* Hartree-Fock energy *) @@ -72,32 +102,28 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64) in (* Convergence criterion *) - let commutator = - let fps = - gemm m_F (gemm m_P m_S) - and spf = - gemm m_S (gemm m_P m_F) - in - Mat.sub fps spf - |> Mat.as_vec + let error = + xt_o_x error_fock m_C + |> Mat.as_vec |> amax + |> abs_float in let converged = - nSCF = max_scf || (abs_float commutator) < threshold_SCF + nSCF = max_scf || error < threshold_SCF in let gap = eigenvalues.{nocc+1} -. eigenvalues.{nocc}; in - Printf.printf "%d %16.10f %11.4e %10.4f\n%!" nSCF energy commutator gap; + Printf.printf "%d %16.10f %11.4e %10.4f\n%!" nSCF energy error gap; if not converged then - loop (nSCF+1) ( (energy, commutator, gap) :: iterations) m_C + loop (nSCF+1) ( (energy, error, gap) :: iterations) m_C diis else let iterations = - List.rev ( (energy, commutator, gap) :: iterations ) + List.rev ( (energy, error, gap) :: iterations ) |> Array.of_list in { HartreeFock_type. @@ -120,6 +146,7 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64) 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 @@ -131,7 +158,8 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64) gemm m_X m_C' in - loop 1 [] m_C + let diis = DIIS.make () in + loop 1 [] m_C diis diff --git a/Simulation.ml b/Simulation.ml index 8766749..f404774 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -12,7 +12,7 @@ type t = { cartesian : bool; } -let make ?cartesian:(cartesian=true) +let make ?cartesian:(cartesian=false) ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei diff --git a/Utils/DIIS.ml b/Utils/DIIS.ml new file mode 100644 index 0000000..fe5c8ba --- /dev/null +++ b/Utils/DIIS.ml @@ -0,0 +1,60 @@ +open Lacaml.D +open Util + +type t = +{ + p : Vec.t list; + e : Vec.t list; + m : int; + mmax : int; +} + +let make ?mmax:(mmax=15) () = + assert (mmax > 1); + { + p = []; + e = []; + m = 0 ; + mmax; + } + + +let append ~p ~e diis = + let update v l = + if diis.m < diis.mmax then + v :: l + else + match List.rev l with + | [] -> assert false + | _ :: rest -> v :: List.rev rest + in + { diis with + p = update p diis.p; + e = update e diis.e; + m = min diis.mmax (diis.m+1); + } + + +let next diis = + let e = Mat.of_col_vecs_list diis.e + and p = Mat.of_col_vecs_list diis.p + in + let a = + let rec aux m = + let a = Mat.make (m+1) (m+1) 1. in + a.{m+1,m+1} <- 0.; + ignore @@ lacpy ~b:a (gemm ~transa:`T ~m ~n:m e e); + if sycon (lacpy a) > 1.e-10 then a + else aux (m-1) + in + aux diis.m + in + let m = Mat.dim1 a - 1 in + let c = Mat.make0 (m+1) 1 in + c.{m+1,1} <- 1.; + sysv a c; + gemm p c ~k:m + |> Mat.as_vec + + + diff --git a/Utils/DIIS.mli b/Utils/DIIS.mli index d3e09c9..bdf057c 100644 --- a/Utils/DIIS.mli +++ b/Utils/DIIS.mli @@ -29,13 +29,13 @@ Equating zero to the derivatives of {% $\mathcal{L}$ %} with respect to {% $c_i$ {% \begin{equation*} \begin{bmatrix} -B_{11} & B_{12} & B_{13} & ... & B_{1m} & -1 \\ -B_{21} & B_{22} & B_{23} & ... & B_{2m} & -1 \\ -B_{31} & B_{32} & B_{33} & ... & B_{3m} & -1 \\ +B_{11} & B_{12} & B_{13} & ... & B_{1m} & 1 \\ +B_{21} & B_{22} & B_{23} & ... & B_{2m} & 1 \\ +B_{31} & B_{32} & B_{33} & ... & B_{3m} & 1 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ -B_{m1} & B_{m2} & B_{m3} & ... & B_{mm} & -1 \\ +B_{m1} & B_{m2} & B_{m3} & ... & B_{mm} & 1 \\ 1 & 1 & 1 & ... & 1 & 0 -\end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_m \\ \lambda \end{bmatrix}= +\end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_m \\ -\lambda \end{bmatrix}= \begin{bmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix} \end{equation*} @@ -50,10 +50,10 @@ $$ %} type t -val make : unit -> t -(** Initialize DIIS *) +val make : ?mmax:int -> unit -> t +(** Initialize DIIS with a maximum size.*) -val append : p:Lacaml.D.Vec.t -> e:Lacaml.D.Vec.t -> t +val append : p:Lacaml.D.Vec.t -> e:Lacaml.D.Vec.t -> t -> t (** Append a parameter vector [p] and the corresponding error vector [e]. *) val next : t -> Lacaml.D.Vec.t diff --git a/Utils/Util.ml b/Utils/Util.ml index a6c186b..be4f8fd 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -171,9 +171,8 @@ let array_product a = let diagonalize_symm m_H = let m_V = lacpy m_H in - let m_W = Vec.create (Mat.dim1 m_H) in let result = - syevd ~vectors:true ~w:m_W m_V + syevd ~vectors:true m_V in m_V, result @@ -200,8 +199,32 @@ let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = gemm c u +let string_of_matrix m = + let open Lacaml.Io in + let rows = Mat.dim1 m + and cols = Mat.dim2 m + in + let rec aux accu first last = + + if (first > last) then String.concat "\n" (List.rev accu) + else + let nw = + Format.asprintf "\n\n %a\n" (Lacaml.Io.pp_lfmat + ~row_labels: + (Array.init rows (fun i -> Printf.sprintf "%d " (i + 1))) + ~col_labels: + (Array.init (min 5 (cols-first+1)) (fun i -> Printf.sprintf "-- %d --" (i + first) )) + ~print_right:false + ~print_foot:false + () ) (lacpy ~ac:first ~n:(min 5 (cols-first+1)) m) + in + aux (nw :: accu) (first+5) last + in + aux [] 1 cols + + let debug_matrix name a = - Format.printf "@[<2>%s =@\n@\n@[%a@]@]@\n@\n" name pp_mat a + Printf.printf "%s =\n%s\n" name (string_of_matrix a)