diff --git a/SCF/RHF.ml b/SCF/RHF.ml index f5ade74..56e1fa4 100644 --- a/SCF/RHF.ml +++ b/SCF/RHF.ml @@ -7,6 +7,31 @@ module El = Electrons module Ao = AOBasis module Ov = Overlap + +type hartree_fock_data = + { + iteration : int ; + coefficients : Mat.t option ; + eigenvalues : Vec.t option ; + error : float option ; + diis : DIIS.t option ; + energy : float option ; + density : Mat.t option ; + fock : Fock.t option ; + } + +let empty = + { + iteration = 0 ; + coefficients = None ; + eigenvalues = None ; + error = None ; + diis = None ; + energy = None ; + density = None ; + fock = None ; + } + let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation = (* Number of occupied MOs *) let nocc = @@ -16,7 +41,7 @@ let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation = let nuclear_repulsion = Si.nuclear_repulsion simulation in - + let ao_basis = Si.ao_basis simulation in @@ -45,138 +70,162 @@ let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation = (* Level shift in MO basis *) let m_LSmo = Array.init (Mat.dim2 m_X) (fun i -> - if i > nocc then level_shift else 0.) + if i > nocc then level_shift else 0.) |> Vec.of_array |> Mat.of_diag in - (* SCF iterations *) - let rec loop nSCF iterations energy_prev m_C m_P_prev fock_prev threshold diis = - (* Density matrix over nocc occupied MOs *) - let m_P = - gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C + (* A single SCF iteration *) + let scf_iteration data = + + let nSCF = data.iteration + 1 + and m_C = of_some data.coefficients + and m_P_prev = data.density + and fock_prev = data.fock + and diis = + match data.diis with + | Some diis -> diis + | None -> DIIS.make () + and threshold = + match data.error with + | Some error -> error + | None -> threshold_SCF *. 2. + in + + (* Density matrix over nocc occupied MOs *) + let m_P = + gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C + in + + (* Fock matrix in AO basis *) + let fock = + match fock_prev, m_P_prev, threshold > 100. *. threshold_SCF with + | Some fock_prev, Some m_P_prev, true -> + let threshold = 1.e-8 in + Fock.make_rhf ~density:(Mat.sub m_P m_P_prev) ~threshold ao_basis + |> Fock.add fock_prev + | _ -> Fock.make_rhf ~density:m_P ao_basis + in + + let m_F, m_Hc, m_J, m_K = + let x = fock in + Fock.(fock x, core x, coulomb x, exchange x) + in + + (* Add level shift in AO basis *) + let m_F = + let m_SC = + gemm m_S m_C in + gemm m_SC (gemm m_LSmo m_SC ~transb:`T) + |> Mat.add m_F + in - (* Fock matrix in AO basis *) - let fock = - match fock_prev, threshold > 100. *. threshold_SCF with - | Some fock_prev, true -> - let threshold = 1.e-8 in - Fock.make_rhf ~density:(Mat.sub m_P m_P_prev) ~threshold ao_basis - |> Fock.add fock_prev - | _ -> Fock.make_rhf ~density:m_P ao_basis + + (* Fock matrix in orthogonal basis *) + let m_F_ortho = + xt_o_x m_F m_X + 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 + xt_o_x (Mat.sub fps spf) m_X + in - let m_F, m_Hc, m_J, m_K = - let x = fock in - Fock.(fock x, core x, coulomb x, exchange x) - in - - (* Add level shift in AO basis *) - let m_F = - let m_SC = - gemm m_S m_C - in - gemm m_SC (gemm m_LSmo m_SC ~transb:`T) - |> Mat.add m_F + let diis = + DIIS.append ~p:(Mat.as_vec m_F_ortho) ~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_ortho) (Mat.dim2 m_F_ortho) + in - (* Fock matrix in orthogonal basis *) - let m_F_ortho = - xt_o_x m_F m_X - in + (* MOs in orthogonal MO basis *) + let m_C', eigenvalues = + diagonalize_symm m_F_diis + 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 - xt_o_x (Mat.sub fps spf) m_X - in + (* MOs in AO basis *) + let m_C = + gemm m_X m_C' + in - let diis = - DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis - in + (* Hartree-Fock energy *) + let energy = + nuclear_repulsion +. 0.5 *. + Mat.gemm_trace m_P (Mat.add m_Hc m_F) + in - let m_F_diis = - let x = - Bigarray.genarray_of_array1 (DIIS.next diis) - in - Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho) - in + (* Convergence criterion *) + let error = + error_fock + |> Mat.as_vec + |> amax + |> abs_float + in + { + iteration = nSCF ; + eigenvalues = Some eigenvalues ; + coefficients = Some m_C ; + error = Some error ; + diis = Some diis ; + energy = Some energy ; + density = Some m_P ; + fock = Some fock ; + } - - (* MOs in orthogonal MO basis *) - let m_C', eigenvalues = - diagonalize_symm m_F_diis - in - - (* MOs in AO basis *) - let m_C = - gemm m_X m_C' - in - - (* Hartree-Fock energy *) - let energy = - nuclear_repulsion +. 0.5 *. - Mat.gemm_trace m_P (Mat.add m_Hc m_F) - in - - (* Convergence criterion *) - let error = - error_fock - |> Mat.as_vec - |> amax - |> abs_float - in - - let converged = - nSCF = max_scf || error < threshold_SCF - in - - let gap = - if nocc < Vec.dim eigenvalues then - eigenvalues.{nocc+1} -. eigenvalues.{nocc} - else 0. - in - - let () = - match energy_prev with - | Some energy_prev -> - Printf.eprintf "%3d %16.10f %16.10f %11.4e %10.4f\n%!" nSCF energy (energy -. energy_prev) error gap - | None -> - Printf.eprintf "%3d %16.10f %16s %11.4e %10.4f\n%!" nSCF energy "" error gap - in - - if not converged then - loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C m_P (Some fock) error diis - else - let iterations = - List.rev ( (energy, error, gap) :: iterations ) - |> Array.of_list - in - HartreeFock_type.(RHF - { - simulation; - nocc; - guess ; - eigenvectors = m_C ; - eigenvalues ; - energy ; - nuclear_repulsion; - iterations ; - kin_energy = Mat.gemm_trace m_P m_T; - 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; - occupation = Mat.copy_diag m_P; - }) in + let rec make_iterations_list data = + + let energy_prev = data.energy in + + (** Perform SCF iteration *) + let data = scf_iteration data in + + (** Check convergence *) + let converged, error = + match data.error with + | None -> false, 0. + | Some error -> (data.iteration = max_scf || error < threshold_SCF), error + in + + (** Print values *) + let nSCF = data.iteration in + + let energy = of_some data.energy in + + let () = + match energy_prev with + | Some energy_prev -> + Printf.eprintf "%3d %16.10f %16.10f %11.4e\n%!" nSCF energy (energy -. energy_prev) error + | None -> + Printf.eprintf "%3d %16.10f %16s %11.4e\n%!" nSCF energy "" error + in + + if converged then + [ data ] + else + { empty with + iteration = data.iteration; + energy = data.energy ; + eigenvalues = data.eigenvalues ; + error = data.error ; + } :: (make_iterations_list data) + in + + + (* Guess coefficients *) let m_H = match guess with @@ -193,8 +242,50 @@ let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation = gemm m_X m_C' in - let diis = DIIS.make () in - loop 1 [] None m_C m_C None threshold_SCF diis + let iterations_list = + make_iterations_list { empty with coefficients = Some m_C } + in + + let iterations, data = + List.map (fun data -> + let gap = + let eigenvalues = of_some data.eigenvalues in + if nocc < Vec.dim eigenvalues then + eigenvalues.{nocc+1} -. eigenvalues.{nocc} + else 0. + and energy = of_some data.energy + and error = of_some data.error + in + (energy, error, gap) + ) iterations_list + |> Array.of_list, + List.hd (List.rev iterations_list) + in + + + let energy = of_some data.energy in + let m_P = of_some data.density in + let fock = of_some data.fock in + let m_J = Fock.coulomb fock in + let m_K = Fock.exchange fock in + + HartreeFock_type.( + RHF { + simulation; + nocc; + guess ; + eigenvectors = of_some data.coefficients ; + eigenvalues = of_some data.eigenvalues ; + energy ; + nuclear_repulsion; + iterations ; + kin_energy = Mat.gemm_trace m_P m_T; + 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; + occupation = Mat.copy_diag m_P; + } + ) diff --git a/Utils/Util.ml b/Utils/Util.ml index 6c7f2f9..ddfd69c 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -171,6 +171,11 @@ let boys_function ~maxm t = end +let of_some = function + | Some a -> a + | None -> assert false + + (** {2 List functions} *) let list_some l = diff --git a/Utils/Util.mli b/Utils/Util.mli index 47c59eb..4971412 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -36,6 +36,7 @@ val chop : float -> (unit -> float) -> float than {!Constants.epsilon}, and return [a *. f ()]. *) +val of_some : 'a option -> 'a (** {2 Functions related to the Boys function} *)