LCPQ
/
QCaml
mirror of https://gitlab.com/scemama/QCaml.git
10
1
Fork 0
Code Releases Activity

Pretty print RHF

This commit is contained in:
Anthony Scemama 2018-05-30 18:07:05 +02:00
parent 25b3c41c50
commit eb54ee8ab5
7 changed files with 122 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
HartreeFock/Fock.ml
View File

@ -4,7 +4,13 @@ open Constants
type t = Mat.t
type t =
{
fock : Mat.t ;
core : Mat.t ;
coulomb : Mat.t ;
exchange : Mat.t ;
}
let make ~density simulation =
let m_P = density
@ -15,24 +21,40 @@ let make ~density simulation =
let nBas = Mat.dim1 m_T
in
let m_F = Mat.add m_T m_V in
let m_Hc = Mat.add m_T m_V
and m_J = Mat.make0 nBas nBas
and m_K = Mat.make0 nBas nBas
in
for sigma = 1 to nBas do
for nu = 1 to nBas do
for lambda = 1 to nBas do
let p = m_P.{lambda,sigma} in
if abs_float p > epsilon then
for mu = 1 to nu do
m_F.{mu,nu} <- m_F.{mu,nu} +. p *.
(ERI.get_phys m_G mu lambda nu sigma
-. 0.5 *. ERI.get_phys m_G mu lambda sigma nu)
m_J.{mu,nu} <- m_J.{mu,nu} +. p *.
ERI.get_phys m_G mu lambda nu sigma
done
done
done
done;
for nu = 1 to nBas do
for sigma = 1 to nBas do
for lambda = 1 to nBas do
let p = 0.5 *. m_P.{lambda,sigma} in
if abs_float p > epsilon then
for mu = 1 to nu do
m_K.{mu,nu} <- m_K.{mu,nu} -. p *.
ERI.get_phys m_G mu lambda sigma nu
done
done
done
done;
for nu = 1 to nBas do
for mu = 1 to nu do
m_F.{nu,mu} <- m_F.{mu,nu}
m_J.{nu,mu} <- m_J.{mu,nu};
m_K.{nu,mu} <- m_K.{mu,nu};
done
done;
m_F
{ fock = Mat.add m_Hc @@ Mat.add m_J m_K ;
core = m_Hc ; coulomb = m_J ; exchange = m_K }

2
HartreeFock/HartreeFock.ml
View File

@ -7,3 +7,5 @@ let make ?guess simulation =
else
invalid_arg "UHF or ROHF not implemented"
let to_string = HartreeFock_type.to_string

43
HartreeFock/HartreeFock_type.ml
View File

@ -1,11 +1,44 @@
type t =
{
guess : Guess.t;
eigenvectors : Lacaml.D.Mat.t ;
eigenvalues : Lacaml.D.Vec.t ;
energy : float ;
iterations : float array;
guess : Guess.t;
eigenvectors : Lacaml.D.Mat.t ;
eigenvalues : Lacaml.D.Vec.t ;
nocc : int;
energy : float ;
nuclear_repulsion : float ;
kin_energy : float ;
eN_energy : float ;
coulomb_energy : float ;
exchange_energy : float ;
iterations : (float * float * float) array;
}
let to_string hf_calc =
let sprintf = Printf.sprintf in
"
=====================================================
Hartree-Fock
=====================================================
# HF energy Convergence HOMO-LUMO
---------------------------------------------------" ::
Array.to_list (Array.mapi (fun i (e, c, g) ->
sprintf " %5d %13.8f %11.4e %11.4f " (i+1) e c g) hf_calc.iterations)
@ [
" =====================================================";
sprintf "%20s energy %16.10f" "One-electron" (hf_calc.kin_energy +. hf_calc.eN_energy) ;
sprintf "%20s energy %16.10f" "Kinetic " hf_calc.kin_energy ;
sprintf "%20s energy %16.10f" "Potential" hf_calc.eN_energy ;
" ---------------------------------------------------";
sprintf "%20s energy %16.10f" "Two-electron" (hf_calc.coulomb_energy +. hf_calc.exchange_energy) ;
sprintf "%20s energy %16.10f" "Coulomb" hf_calc.coulomb_energy ;
sprintf "%20s energy %16.10f" "Exchange" hf_calc.exchange_energy ;
" ---------------------------------------------------";
sprintf "%20s energy %16.10f" "Electronic" (hf_calc.energy -. hf_calc.nuclear_repulsion);
sprintf "%17s repulsion %16.10f" "Nuclear" hf_calc.coulomb_energy ;
sprintf "%20s energy %16.10f" "Hartree-Fock" hf_calc.energy ;
" =====================================================" ]
|> String.concat "\n"

19
HartreeFock/HartreeFock_type.mli
View File

@ -2,11 +2,20 @@
type t =
{
guess : Guess.t; (** Initial guess *)
eigenvectors : Lacaml.D.Mat.t ; (** Final eigenvectors *)
eigenvalues : Lacaml.D.Vec.t ; (** Final eigenvalues *)
energy : float ; (** Final energy *)
iterations : float array; (** Energies of all iterations *)
guess : Guess.t; (** Initial guess *)
eigenvectors : Lacaml.D.Mat.t ; (** Final eigenvectors *)
eigenvalues : Lacaml.D.Vec.t ; (** Final eigenvalues *)
nocc : int ; (** Number of occupied MOs *)
energy : float ; (** Final energy *)
nuclear_repulsion : float ; (** Nucleus-Nucleus potential energy *)
kin_energy : float ; (** Kinetic energy *)
eN_energy : float ; (** Electron-nucleus potential energy *)
coulomb_energy : float ; (** Electron-Electron potential energy *)
exchange_energy : float ; (** Exchange energy *)
iterations : (float * float * float) array;
(** Energy, convergence and HOMO-LUMO gap of all iterations *)
}
val to_string : t -> string
(** Results of a Hartree-Fock calculation pretty-printed in a string. *)

43
HartreeFock/RHF.ml
View File

@ -10,6 +10,10 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
simulation.electrons.Electrons.n_alpha
in
let nuclear_repulsion =
simulation.nuclear_repulsion
in
(* Initial guess *)
let guess =
Guess.make ~guess simulation
@ -21,18 +25,14 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
in
(* Core Hamiltonian *)
let m_Hc =
let m_T = Lazy.force simulation.kin_ints
and m_V = Lazy.force simulation.eN_ints
in Mat.add m_T m_V
in
(* Overlap matrix *)
let m_S =
Lazy.force simulation.overlap
in
let m_T = Lazy.force simulation.kin_ints
and m_V = Lazy.force simulation.eN_ints
in
(* SCF iterations *)
let rec loop nSCF iterations m_C =
@ -43,8 +43,11 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
in
(* Fock matrix in AO basis *)
let m_F =
Fock.make ~density:m_P simulation
let m_F, m_Hc, m_J, m_K =
let x =
Fock.make ~density:m_P simulation
in
x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange
in
(* Fock matrix in MO basis *)
@ -64,7 +67,7 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
(* Hartree-Fock energy *)
let energy =
simulation.nuclear_repulsion +. 0.5 *.
nuclear_repulsion +. 0.5 *.
Mat.gemm_trace m_P (Mat.add m_Hc m_F)
in
@ -81,24 +84,34 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
in
let converged =
nSCF = max_scf || (abs_float commutator) < threshold_SCF
nSCF = max_scf || (abs_float commutator) < threshold_SCF
in
Printf.printf "%d %16.10f %10.4e\n%!" nSCF energy commutator;
let gap =
eigenvalues.{nocc+1} -. eigenvalues.{nocc};
in
Printf.printf "%d %16.10f %11.4e %10.4f\n%!" nSCF energy commutator gap;
if not converged then
loop (nSCF+1) (energy :: iterations) m_C
loop (nSCF+1) ( (energy, commutator, gap) :: iterations) m_C
else
let iterations =
List.rev iterations
List.rev ( (energy, commutator, gap) :: iterations )
|> Array.of_list
in
{ HartreeFock_type.
nocc;
guess ;
eigenvectors = m_C ;
eigenvalues ;
energy = iterations.(Array.length iterations - 1) ;
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;
}
in

12
Utils/Util.mli
View File

@ -58,7 +58,10 @@ val array_product : float array -> float
(** {2 Extension of the List module} *)
val list_some : 'a option list -> 'a list
(** Filters out all [None] elements of the list, and returns the elements without
the [Some]. *)
(** {2 Linear algebra } *)
@ -72,11 +75,12 @@ val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat
val canonical_ortho: ?thresh:float -> overlap:Lacaml.D.mat -> Lacaml.D.mat -> Lacaml.D.mat
(** Canonical orthogonalization. [overlap] is the overlap matrix {% $\mathbf{S}$ %},
and the last argument contains the vectors {% $\mathbf{C}$ %} to orthogonalize.
{%
\begin{align}
\mathbf{U\, D\, V^\dag} & = \mathbf{S} \\
\mathbf{C_\bot} & = \mathbf{C\, U\, D^{-1/2}}
\end{align}
$$
\mathbf{C_\bot} = \mathbf{C\, U\, D^{-1/2}}, \;
\mathbf{U\, D\, V^\dag} = \mathbf{S}
$$
%}
*)

8
run_hartree_fock.ml
View File

@ -39,11 +39,9 @@ let run ~out =
Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file
in
let _ =
HartreeFock.make s
in
Printf.printf "Done.\n%!";
()
HartreeFock.make s
|> HartreeFock.to_string
|> print_endline
let () =