LCPQ/QCaml
10
1
Fork 0
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-12-26 22:33:36 +01:00
Code Releases Activity

Printing of HF MOs

Browse Source
This commit is contained in:
Anthony Scemama 2018-05-30 20:16:18 +02:00
parent eb54ee8ab5
commit 07ee916a28
3 changed files with 103 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

116
HartreeFock/HartreeFock_type.ml
View File

@ -1,3 +1,6 @@
open Lacaml.D
open Util
type t =
{
guess : Guess.t;
@ -14,31 +17,98 @@ type t =
}
let to_string hf_calc =
let sprintf = Printf.sprintf in
"
=====================================================
Hartree-Fock
=====================================================
# HF energy Convergence HOMO-LUMO
let iterations_to_string hf_calc =
" # 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 ;
" =====================================================" ]
Printf.sprintf " %5d %13.8f %11.4e %11.4f " (i+1) e c g) hf_calc.iterations)
@ [ " -----------------------------------------------------" ]
|> String.concat "\n"
let summary hf_calc =
let ha_to_ev = Constants.ha_to_ev in
String.concat "\n" [
" =====================================================";
Printf.sprintf " One-electron energy %16.10f" (hf_calc.kin_energy +. hf_calc.eN_energy) ;
Printf.sprintf " Kinetic energy %16.10f" hf_calc.kin_energy ;
Printf.sprintf " Potential energy %16.10f" hf_calc.eN_energy ;
" ---------------------------------------------------";
Printf.sprintf " Two-electron energy %16.10f" (hf_calc.coulomb_energy +. hf_calc.exchange_energy) ;
Printf.sprintf " Coulomb energy %16.10f" hf_calc.coulomb_energy ;
Printf.sprintf " Exchange energy %16.10f" hf_calc.exchange_energy ;
" ---------------------------------------------------";
Printf.sprintf " Electronic energy %16.10f" (hf_calc.energy -. hf_calc.nuclear_repulsion);
Printf.sprintf " Nuclear repulsion %16.10f" hf_calc.coulomb_energy ;
Printf.sprintf " Hartree-Fock energy %16.10f" hf_calc.energy ;
" ---------------------------------------------------";
Printf.sprintf " HF HOMO energy %16.10f eV" (ha_to_ev *. hf_calc.eigenvalues.{hf_calc.nocc} );
Printf.sprintf " HF LUMO energy %16.10f eV" (ha_to_ev *. hf_calc.eigenvalues.{hf_calc.nocc+1} );
Printf.sprintf " HF HOMO-LUMO gap %16.10f eV"
(ha_to_ev *.(hf_calc.eigenvalues.{hf_calc.nocc+1} -. hf_calc.eigenvalues.{hf_calc.nocc}));
" =====================================================" ]
let mos_to_string hf_calc =
let open Lacaml.Io in
let rows = Mat.dim1 hf_calc.eigenvectors
and cols = Mat.dim2 hf_calc.eigenvectors
in
let rec aux accu first last =
if (first > last) then String.concat "\n" (List.rev accu)
else
let nw =
"\n Eigenvalues : " ^ (
Array.mapi (fun i x ->
if (i+1 >= first) && (i+1 <= first+4 ) then
Some (Printf.sprintf " %f " x)
else None
)
(Vec.to_array hf_calc.eigenvalues)
|> Array.to_list
|> list_some
|> String.concat " "
)^
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)) hf_calc.eigenvectors)
in
aux (nw :: accu) (first+5) last
in
String.concat "\n" [ "
=========================================================================
Molecular Orbitals
=========================================================================
Occupied
-----------------------------------------------------------------------
" ;
aux [] 1 hf_calc.nocc ;
"
Virtual
-----------------------------------------------------------------------
" ;
aux [] (hf_calc.nocc+1) (Mat.dim2 hf_calc.eigenvectors) ;
"
=========================================================================
" ]
let to_string hf_calc =
String.concat "\n" [ "
=====================================================
Hartree-Fock
=====================================================" ; "" ;
iterations_to_string hf_calc ; "" ;
summary hf_calc ; "" ;
mos_to_string hf_calc ; "" ;
]

4
Utils/Constants.ml
View File

@ -8,5 +8,7 @@ let sq_pi_over_two = sq_pi *. 0.5
let pi_inv = 1. /. pi
let two_over_sq_pi = 2. /. (sqrt pi)
let a0 = 0.529_177_210_671_2
let a0 = 0.529_177_210_67
let a0_inv = 1. /. a0
let ha_to_ev = 27.211_386_02
let ev_to_ha = 1. /. ha_to_ev

8
Utils/Constants.mli
View File

@ -29,8 +29,14 @@ val two_over_sq_pi : float
(** {2 Physical constants} *)
val a0 : float
(** Bohr radius : {% $a_0$ %} = 0.529_177_210_671_2 Angstrom *)
(** Bohr radius : {% $a_0$ %} = 0.529 177 210 67(23) Angstrom *)
val a0_inv : float
(** {% $\frac{1}{a_0}$ %} *)
val ha_to_ev : float
(** Hartree to eV conversion factor :: 27.211 386 02(17) *)
val ev_to_ha : float
(** eV to Hartree conversion factor : 1/{ha_to_ev} *)