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

A lot of cleaning in HF

This commit is contained in:
Anthony Scemama 2019-03-01 21:56:46 +01:00
parent 4b9e2016aa
commit 6e68eef645
11 changed files with 837 additions and 780 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
MOBasis/MOBasis.ml
View File

@ -4,11 +4,11 @@ open Constants
(** One-electron orthogonal basis set, corresponding to Molecular Orbitals. *)
module HF = HartreeFock_type
module HF = HartreeFock
module Si = Simulation
type mo_type =
| RHF | ROHF | CASSCF
| RHF | ROHF | UHF | CASSCF
| Natural of string
| Localized of string
@ -39,6 +39,22 @@ let kin_ints t = Lazy.force t.kin_ints
let two_e_ints t = Lazy.force t.ee_ints
let one_e_ints t = Lazy.force t.one_e_ints
let mo_energies t =
let f =
let m_C = mo_coef t in
let m_N = Mat.of_diag @@ mo_occupation t in
let m_P =
gemm m_C @@ (gemm m_N m_C ~transb:`T)
in
match t.mo_type with
| RHF -> Fock.make_rhf ~density:m_P (ao_basis t)
| ROHF -> Fock.make_uhf ~density_same:m_P ~density_other:m_P (ao_basis t)
| _ -> failwith "Not implemented"
in
let m_F = Fock.fock f in
Vec.init (size t) (fun i -> m_F.{i,i})
let mo_matrix_of_ao_matrix ~mo_coef ao_matrix =
xt_o_x ~x:mo_coef ~o:ao_matrix
@ -54,8 +70,8 @@ let four_index_transform ~mo_coef eri_ao =
let mo_num = Mat.dim2 mo_coef in
let eri_mo = ERI.create ~size:mo_num `Dense in
let mo_num_2 = mo_num * mo_num in
let ao_num_2 = ao_num * ao_num in
let mo_num_2 = mo_num * mo_num in
let ao_num_2 = ao_num * ao_num in
let ao_mo_num = ao_num * mo_num in
let range_mo = list_range 1 mo_num in
@ -156,26 +172,58 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef () =
eN_ints ; ee_ints ; kin_ints ; one_e_ints }
let of_rhf hf =
let mo_coef = hf.HF.eigenvectors in
let simulation = hf.HF.simulation in
let mo_type = RHF in
let mo_occupation = hf.HF.occupation in
let of_hartree_fock hf =
let mo_coef = HF.eigenvectors hf in
let simulation = HF.simulation hf in
let mo_occupation = HF.occupation hf in
let mo_type =
match HF.kind hf with
| HartreeFock.RHF -> RHF
| HartreeFock.ROHF -> ROHF
| HartreeFock.UHF -> UHF
in
make ~simulation ~mo_type ~mo_occupation ~mo_coef ()
let of_rohf hf =
let mo_coef = hf.HF.eigenvectors in
let simulation = hf.HF.simulation in
let mo_type = ROHF in
let mo_occupation = hf.HF.occupation in
make ~simulation ~mo_type ~mo_occupation ~mo_coef ()
let of_hartree_fock = function
| HF.RHF hf -> of_rhf hf
| HF.ROHF hf -> of_rohf hf
| HF.UHF _ -> assert false
let pp_mo ?(start=1) ?finish ppf t =
let open Lacaml.Io in
let rows = Mat.dim1 t.mo_coef
and cols = Mat.dim2 t.mo_coef
in
let finish =
match finish with
| None -> cols
| Some x -> x
in
let rec aux first =
if (first > finish) then ()
else
begin
Format.fprintf ppf "@[<v>@[<v4>@[<h>%s@;" "Eigenvalues:";
Array.iteri (fun i x ->
if (i+1 >= first) && (i+1 <= first+4 ) then
Format.fprintf ppf "%12f@ " x)
(Vec.to_array @@ mo_energies t);
Format.fprintf ppf "@]@;";
Format.fprintf ppf "@[%a@]"
(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)) (t.mo_coef)) ;
Format.fprintf ppf "@]@;@;@]";
aux (first+5)
end
in
aux start

9
MOBasis/MOBasis.mli
View File

@ -7,7 +7,7 @@
open Lacaml.D
type mo_type =
| RHF | ROHF | CASSCF
| RHF | ROHF | UHF | CASSCF
| Natural of string
| Localized of string
@ -58,7 +58,12 @@ val make : simulation:Simulation.t ->
unit -> t
(** Function to build a data structure representing the molecular orbitals. *)
val of_hartree_fock : HartreeFock_type.t -> t
val of_hartree_fock : HartreeFock.t -> t
(** Build MOs from a Restricted Hartree-Fock calculation. *)
(** {1 Printers} *)
val pp_mo : ?start:int -> ?finish:int -> Format.formatter -> t -> unit

678
SCF/HartreeFock.ml
View File

@ -1,18 +1,678 @@
open Lacaml.D
open Util
open Constants
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 ;
density_a : Mat.t option ;
density_b : Mat.t option ;
fock : Fock.t option ;
fock_a : Fock.t option ;
fock_b : Fock.t option ;
}
type hartree_fock_kind =
| RHF (** Restricted Hartree-Fock *)
| ROHF (** Restricted Open-shell Hartree-Fock *)
| UHF (** Unrestricted Hartree-Fock *)
type t =
{
kind : hartree_fock_kind;
simulation : Simulation.t;
guess : Guess.t;
data : hartree_fock_data option lazy_t array;
nocc : int ;
}
let empty =
{
iteration = 0 ;
coefficients = None ;
eigenvalues = None ;
error = None ;
diis = None ;
energy = None ;
density = None ;
density_a = None ;
density_b = None ;
fock = None ;
fock_a = None ;
fock_b = None ;
}
module Si = Simulation
module El = Electrons
module Ao = AOBasis
module Ov = Overlap
let kind t = t.kind
let simulation t = t.simulation
let guess t = t.guess
let nocc t = t.nocc
let n_iterations t =
Array.fold_left (fun accu x ->
match Lazy.force x with
| Some x -> accu + 1
| None -> accu
) 0 t.data
let last_iteration t =
of_some @@ Lazy.force (t.data.(n_iterations t - 1))
let eigenvectors t =
let data = last_iteration t in
of_some data.coefficients
let eigenvalues t =
let data = last_iteration t in
of_some data.eigenvalues
let density t =
let data = last_iteration t in
match kind t with
| RHF -> of_some data.density
| ROHF -> Mat.add (of_some data.density_a) (of_some data.density_b)
| _ -> failwith "Not implemented"
let occupation t =
density t
|> Mat.copy_diag ~n:(Mat.dim2 @@ eigenvectors t)
let energy t =
let data = last_iteration t in
of_some data.energy
let nuclear_repulsion t =
Si.nuclear_repulsion (simulation t)
let ao_basis t =
Si.ao_basis (simulation t)
let kin_energy t =
let m_T =
ao_basis t
|> Ao.kin_ints
|> KinInt.matrix
in
let m_P = density t in
Mat.gemm_trace m_P m_T
let eN_energy t =
let m_V =
ao_basis t
|> Ao.eN_ints
|> NucInt.matrix
in
let m_P = density t in
Mat.gemm_trace m_P m_V
let coulomb_energy t =
let data =
last_iteration t
in
match kind t with
| RHF -> let m_P = of_some data.density in
let fock = of_some data.fock in
let m_J = Fock.coulomb fock in
0.5 *. Mat.gemm_trace m_P m_J
| ROHF -> let m_P_a = of_some data.density_a in
let m_P_b = of_some data.density_b in
let fock_a = of_some data.fock_a in
let fock_b = of_some data.fock_b in
let m_J_a = Fock.coulomb fock_a in
let m_J_b = Fock.coulomb fock_b in
0.5 *. ( (Mat.gemm_trace m_P_a m_J_a) +. (Mat.gemm_trace m_P_b m_J_b) )
| _ -> failwith "Not implemented"
let exchange_energy t =
let data =
last_iteration t
in
match kind t with
| RHF -> let m_P = of_some data.density in
let fock = of_some data.fock in
let m_K = Fock.exchange fock in
0.5 *. Mat.gemm_trace m_P m_K
| ROHF -> let m_P_a = of_some data.density_a in
let m_P_b = of_some data.density_b in
let fock_a = of_some data.fock_a in
let fock_b = of_some data.fock_b in
let m_K_a = Fock.exchange fock_a in
let m_K_b = Fock.exchange fock_b in
0.5 *. ( (Mat.gemm_trace m_P_a m_K_a) +. (Mat.gemm_trace m_P_b m_K_b) )
| _ -> failwith "Not implemented"
let make
?kind
?guess:(guess=`Huckel)
?max_scf:(max_scf=64)
?level_shift:(level_shift=0.1)
?threshold_SCF:(threshold_SCF=1.e-8)
simulation =
let f =
if Electrons.multiplicity @@ Simulation.electrons simulation = 1 then
RHF.make
else
ROHF.make
in f ~guess ~max_scf ~level_shift ~threshold_SCF simulation
simulation =
let to_string = HartreeFock_type.to_string
(* Number of occupied MOs *)
let n_alfa, n_beta =
El.n_alfa @@ Si.electrons simulation,
El.n_beta @@ Si.electrons simulation
in
let nocc = n_alfa in
let kind =
match kind with
| Some kind -> kind
| None -> if (n_alfa = n_beta) then RHF else ROHF
in
let nuclear_repulsion =
Si.nuclear_repulsion simulation
in
let ao_basis =
Si.ao_basis simulation
in
(* Initial guess *)
let guess =
Guess.make ~nocc ~guess ao_basis
in
(* Orthogonalization matrix *)
let m_X =
Ao.ortho ao_basis
in
(* Overlap matrix *)
let m_S =
Ao.overlap ao_basis
|> Ov.matrix
in
(* Level shift in MO basis *)
let m_LSmo =
Array.init (Mat.dim2 m_X) (fun i ->
if i > nocc then level_shift else 0.)
|> Vec.of_array
|> Mat.of_diag
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'
in
(* A single SCF iteration *)
let scf_iteration_rhf 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 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 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
(* 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
{ empty with
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 ;
}
in
let scf_iteration_rohf data =
let nSCF = data.iteration + 1
and m_C = of_some data.coefficients
and m_P_a_prev = data.density_a
and m_P_b_prev = data.density_b
and fock_a_prev = data.fock_a
and fock_b_prev = data.fock_b
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 *)
let m_P_a =
gemm ~alpha:1. ~transb:`T ~k:n_alfa m_C m_C
in
let m_P_b =
gemm ~alpha:1. ~transb:`T ~k:n_beta m_C m_C
in
let m_P =
Mat.add m_P_a m_P_b
in
(* Fock matrix in AO basis *)
let fock_a =
match fock_a_prev, threshold > 100. *. threshold_SCF with
| Some fock_a_prev, true ->
let threshold = 1.e-8 in
Fock.make_uhf ~density_same:(Mat.sub m_P_a @@ of_some m_P_a_prev) ~density_other:(Mat.sub m_P_b @@ of_some m_P_b_prev) ~threshold ao_basis
|> Fock.add fock_a_prev
| _ -> Fock.make_uhf ~density_same:m_P_a ~density_other:m_P_b ao_basis
in
let fock_b =
match fock_b_prev, threshold > 100. *. threshold_SCF with
| Some fock_b_prev, true ->
let threshold = 1.e-8 in
Fock.make_uhf ~density_same:(Mat.sub m_P_b @@ of_some m_P_b_prev) ~density_other:(Mat.sub m_P_a @@ of_some m_P_a_prev) ~threshold ao_basis
|> Fock.add fock_b_prev
| _ -> Fock.make_uhf ~density_same:m_P_b ~density_other:m_P_a ao_basis
in
let m_F_a, m_Hc_a, m_J_a, m_K_a =
let x = fock_a in
Fock.(fock x, core x, coulomb x, exchange x)
in
let m_F_b, m_Hc_b, m_J_b, m_K_b =
let x = fock_b in
Fock.(fock x, core x, coulomb x, exchange x)
in
let m_F_mo_a =
xt_o_x ~o:m_F_a ~x:m_C
in
let m_F_mo_b =
xt_o_x ~o:m_F_b ~x:m_C
in
let m_F_mo =
let space k =
if k <= n_beta then
`Core
else if k <= n_alfa then
`Active
else
`Virtual
in
Array.init (Mat.dim2 m_F_mo_a) (fun i ->
let i = i+1 in
Array.init (Mat.dim1 m_F_mo_a) (fun j ->
let j = j+1 in
match (space i), (space j) with
| `Core , `Core ->
0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) -.
(m_F_mo_b.{i,j} -. m_F_mo_a.{i,j})
| `Active , `Core
| `Core , `Active ->
m_F_mo_b.{i,j}
| `Core , `Virtual
| `Virtual , `Core
| `Active , `Active ->
0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j})
| `Virtual , `Active
| `Active , `Virtual ->
m_F_mo_a.{i,j}
| `Virtual , `Virtual ->
0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) +.
(m_F_mo_b.{i,j} -. m_F_mo_a.{i,j})
) )
|> Mat.of_array
in
let m_SC =
gemm m_S m_C
in
let m_F =
x_o_xt ~x:m_SC ~o:m_F_mo
in
(* Add level shift in AO basis *)
let m_F =
x_o_xt ~x:m_SC ~o:m_LSmo
|> Mat.add m_F
in
(* 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 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
(* 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_a (Mat.add m_Hc_a m_F_a) +.
Mat.gemm_trace m_P_b (Mat.add m_Hc_b m_F_b) )
in
(* Convergence criterion *)
let error =
error_fock
|> Mat.as_vec
|> amax
|> abs_float
in
{ empty with
iteration = nSCF ;
eigenvalues = Some eigenvalues ;
coefficients = Some m_C ;
error = Some error ;
diis = Some diis ;
energy = Some energy ;
density_a = Some m_P_a ;
density_b = Some m_P_b ;
fock_a = Some fock_a ;
fock_b = Some fock_b ;
}
in
let scf_iteration data =
match kind with
| RHF -> scf_iteration_rhf data
| ROHF -> scf_iteration_rohf data
| _ -> failwith "Not implemented"
in
let array_data =
let storage =
Array.make max_scf None
in
let rec iteration = function
| 0 -> Some (scf_iteration { empty with coefficients = Some m_C })
| n -> begin
match storage.(n) with
| Some result -> Some result
| None ->
begin
let data =
(iteration (n-1))
in
match data with
| None -> None
| Some data ->
begin
(** Check convergence *)
let converged, error =
match data.error with
| None -> false, 0.
| Some error -> (data.iteration = max_scf || error < threshold_SCF), error
in
if converged then
None
else
begin
storage.(n-1) <- Some { empty with
iteration = data.iteration;
energy = data.energy ;
eigenvalues = data.eigenvalues ;
error = data.error ;
};
storage.(n) <- Some (scf_iteration data);
storage.(n);
end
end
end
end
in
Array.init max_scf (fun i -> lazy (iteration i))
in
{
kind;
simulation;
guess ;
data = array_data;
nocc;
}
let linewidth = 60
let pp_iterations ppf t =
Format.fprintf ppf "@[%4s%s@]@." "" (Printing.line linewidth);
Format.fprintf ppf "@[%4s@[%5s@]@,@[%16s@]@,@[%16s@]@,@[%11s@]@]@."
"" "#" "HF energy " "Convergence" "HOMO-LUMO";
Format.fprintf ppf "@[%4s%s@]@." "" (Printing.line linewidth);
let nocc = nocc t in
Array.iter (fun data ->
let data = Lazy.force data in
match data with
| None -> ()
| Some data ->
let e = of_some data.eigenvalues in
let gap = e.{nocc+1} -. e.{nocc} in
begin
Format.fprintf ppf "@[%4s@[%5d@]@,@[%16.8f@]@,@[%16.4e@]@,@[%11.4f@]@]@." ""
(data.iteration) (of_some data.energy) (of_some data.error) gap;
end
) t.data;
Format.fprintf ppf "@[%4s%s@]@." "" (Printing.line linewidth)
let pp_summary ppf t =
let print text value =
Format.fprintf ppf "@[@[%30s@]@,@[%16.10f@]@]@;" text value;
and line () =
Format.fprintf ppf "@[ %s @]@;" (Printing.line (linewidth-4));
in
let ha_to_ev = Constants.ha_to_ev in
let e = eigenvalues t in
Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth);
Format.fprintf ppf "@[<v>";
print "One-electron energy" (kin_energy t +. eN_energy t) ;
print "Kinetic" (kin_energy t) ;
print "Potential" (eN_energy t) ;
line () ;
print "Two-electron energy" (coulomb_energy t +. exchange_energy t) ;
print "Coulomb" (coulomb_energy t) ;
print "Exchange" (exchange_energy t) ;
line ();
print "HF HOMO" (ha_to_ev *. e.{nocc t});
print "HF LUMO" (ha_to_ev *. e.{nocc t + 1});
print "HF LUMO-LUMO" (ha_to_ev *. (e.{nocc t + 1} -. e.{nocc t }));
line ();
print "Electronic energy" (energy t -. nuclear_repulsion t) ;
print "Nuclear repulsion" (nuclear_repulsion t) ;
print "Hartree-Fock energy" (energy t) ;
Format.fprintf ppf "@]";
Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth)
let pp_hf ppf t =
Format.fprintf ppf "@.@[%s@]@." (Printing.line ~c:'=' 70);
Format.fprintf ppf "@[%34s %-34s@]@." (match t.kind with
| UHF -> "Unrestricted"
| RHF -> "Restricted"
| ROHF -> "Restricted Open-shell") "Hartree-Fock";
Format.fprintf ppf "@[%s@]@.@." (Printing.line ~c:'=' 70);
Format.fprintf ppf "@[%a@]@." pp_iterations t;
Format.fprintf ppf "@[<v 4>@;%a@]@." pp_summary t

76
SCF/HartreeFock.mli Normal file
View File

@ -0,0 +1,76 @@
open Lacaml.D
(** Data structure representing the output of a Hartree-Fock caculation *)
type hartree_fock_data
type hartree_fock_kind =
| RHF (** Restricted Hartree-Fock *)
| ROHF (** Restricted Open-shell Hartree-Fock *)
| UHF (** Unrestricted Hartree-Fock *)
type t
val kind : t -> hartree_fock_kind
(** Kind of simulation : RHF, ROHF or UHF. *)
val simulation : t -> Simulation.t
(** Simulation which was used for HF calculation *)
val guess : t -> Guess.t
(** Initial guess *)
val eigenvectors : t -> Mat.t
(** Final eigenvectors *)
val eigenvalues : t -> Vec.t
(** Final eigenvalues *)
val occupation : t -> Vec.t
(** Diagonal of the density matrix *)
val energy : t -> float
(** Final energy *)
val nuclear_repulsion : t -> float
(** Nucleus-Nucleus potential energy *)
val kin_energy : t -> float
(** Kinetic energy *)
val eN_energy : t -> float
(** Electron-nucleus potential energy *)
val coulomb_energy : t -> float
(** Electron-Electron potential energy *)
val exchange_energy : t -> float
(** Exchange energy *)
val nocc : t -> int
(** Number of occupied MOs *)
val empty: hartree_fock_data
(** Empty data *)
val make :
?kind:hartree_fock_kind ->
?guess:[ `Hcore | `Huckel ] ->
?max_scf:int ->
?level_shift:float -> ?threshold_SCF:float -> Simulation.t -> t
(** {1 Printers} *)
val pp_hf : Format.formatter -> t -> unit
val pp_iterations : Format.formatter -> t -> unit
val pp_summary : Format.formatter -> t -> unit
(*
val pp_mos : Format.formatter -> t -> unit
#*)

128
SCF/HartreeFock_type.ml
View File

@ -1,128 +0,0 @@
open Lacaml.D
open Util
type s =
{
simulation : Simulation.t;
guess : Guess.t;
eigenvectors : Mat.t ;
eigenvalues : Vec.t ;
occupation : Vec.t ;
iterations : (float * float * float) array;
energy : float ;
nuclear_repulsion : float ;
kin_energy : float ;
eN_energy : float ;
coulomb_energy : float ;
exchange_energy : float ;
nocc : int;
}
type t =
| RHF of s (** Restricted Hartree-Fock *)
| ROHF of s (** Restricted Open-shell Hartree-Fock *)
| UHF of s (** Unrestricted Hartree-Fock *)
let iterations_to_string hf_calc =
" # HF energy Convergence HOMO-LUMO
---------------------------------------------------" ::
Array.to_list (Array.mapi (fun i (e, c, g) ->
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 =
let aux hf_calc r =
String.concat "\n" [ Printf.sprintf "
=====================================================
%s Hartree-Fock
=====================================================" r ; "" ;
iterations_to_string hf_calc ; "" ;
summary hf_calc ; "" ;
mos_to_string hf_calc ; "" ;
]
in
match hf with
| RHF hf_calc -> aux hf_calc "Restricted"
| UHF hf_calc -> aux hf_calc "Unrestricted"
| ROHF hf_calc -> aux hf_calc "Restricted Open-shell"

31
SCF/HartreeFock_type.mli
View File

@ -1,31 +0,0 @@
open Lacaml.D
(** Data structure representing the output of a Hartree-Fock caculation *)
type s =
{
simulation : Simulation.t; (** Simulation which was used for HF calculation *)
guess : Guess.t; (** Initial guess *)
eigenvectors : Mat.t ; (** Final eigenvectors *)
eigenvalues : Vec.t ; (** Final eigenvalues *)
occupation : Vec.t ; (** Diagonal of the density matrix *)
iterations : (float * float * float) array;
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 *)
nocc : int ; (** Number of occupied MOs *)
(** Energy, convergence and HOMO-LUMO gap of all iterations *)
}
type t =
| RHF of s (** Restricted Hartree-Fock *)
| ROHF of s (** Restricted Open-shell Hartree-Fock *)
| UHF of s (** Unrestricted Hartree-Fock *)
val to_string : t -> string
(** Results of a Hartree-Fock calculation pretty-printed in a string. *)

291
SCF/RHF.ml
View File

@ -1,291 +0,0 @@
open Util
open Constants
open Lacaml.D
module Si = Simulation
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 =
El.n_alfa @@ Si.electrons simulation
in
let nuclear_repulsion =
Si.nuclear_repulsion simulation
in
let ao_basis =
Si.ao_basis simulation
in
(* Initial guess *)
let guess =
Guess.make ~nocc ~guess ao_basis
in
(* Orthogonalization matrix *)
let m_X =
Ao.ortho ao_basis
in
(* Overlap matrix *)
let m_S =
Ao.overlap ao_basis
|> Ov.matrix
in
let m_T = Ao.kin_ints ao_basis |> KinInt.matrix
and m_V = Ao.eN_ints ao_basis |> NucInt.matrix
in
(* Level shift in MO basis *)
let m_LSmo =
Array.init (Mat.dim2 m_X) (fun i ->
if i > nocc then level_shift else 0.)
|> Vec.of_array
|> Mat.of_diag
in
(* 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 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 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
(* 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
{
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 ;
}
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
| 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
let m_C =
gemm m_X m_C'
in
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;
}
)

290
SCF/ROHF.ml
View File

@ -1,290 +0,0 @@
open Util
open Constants
open Lacaml.D
module Si = Simulation
module El = Electrons
module Ao = AOBasis
module Ov = Overlap
let make ~guess ~max_scf ~level_shift ~threshold_SCF simulation =
(* Number of occupied MOs *)
let n_alfa =
El.n_alfa @@ Si.electrons simulation
in
let nocc = n_alfa in
let n_beta =
El.n_beta @@ Si.electrons simulation
in
let nuclear_repulsion =
Si.nuclear_repulsion simulation
in
let ao_basis =
Si.ao_basis simulation
in
(* Initial guess *)
let guess =
Guess.make ~nocc ~guess ao_basis
in
(* Orthogonalization matrix *)
let m_X =
Ao.ortho ao_basis
in
(* Overlap matrix *)
let m_S =
Ao.overlap ao_basis
|> Ov.matrix
in
let m_T = Ao.kin_ints ao_basis |> KinInt.matrix
and m_V = Ao.eN_ints ao_basis |> NucInt.matrix
in
(* Level shift in MO basis *)
let m_LSmo =
Array.init (Mat.dim2 m_X) (fun i ->
if i > n_alfa 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_a_prev m_P_b_prev fock_a_prev fock_b_prev threshold diis =
(* Density matrix *)
let m_P_a =
gemm ~alpha:1. ~transb:`T ~k:n_alfa m_C m_C
in
let m_P_b =
gemm ~alpha:1. ~transb:`T ~k:n_beta m_C m_C
in
let m_P =
Mat.add m_P_a m_P_b
in
(* Fock matrix in AO basis *)
let fock_a =
match fock_a_prev, threshold > 100. *. threshold_SCF with
| Some fock_a_prev, true ->
let threshold = 1.e-8 in
Fock.make_uhf ~density_same:(Mat.sub m_P_a m_P_a_prev) ~density_other:(Mat.sub m_P_b m_P_b_prev) ~threshold ao_basis
|> Fock.add fock_a_prev
| _ -> Fock.make_uhf ~density_same:m_P_a ~density_other:m_P_b ao_basis
in
let fock_b =
match fock_b_prev, threshold > 100. *. threshold_SCF with
| Some fock_b_prev, true ->
let threshold = 1.e-8 in
Fock.make_uhf ~density_same:(Mat.sub m_P_b m_P_b_prev) ~density_other:(Mat.sub m_P_a m_P_a_prev) ~threshold ao_basis
|> Fock.add fock_b_prev
| _ -> Fock.make_uhf ~density_same:m_P_b ~density_other:m_P_a ao_basis
in
let m_F_a, m_Hc_a, m_J_a, m_K_a =
let x = fock_a in
Fock.(fock x, core x, coulomb x, exchange x)
in
let m_F_b, m_Hc_b, m_J_b, m_K_b =
let x = fock_b in
Fock.(fock x, core x, coulomb x, exchange x)
in
let m_F_mo_a =
xt_o_x ~o:m_F_a ~x:m_C
in
let m_F_mo_b =
xt_o_x ~o:m_F_b ~x:m_C
in
let m_F_mo =
let space k =
if k <= n_beta then
`Core
else if k <= n_alfa then
`Active
else
`Virtual
in
Array.init (Mat.dim2 m_F_mo_a) (fun i ->
let i = i+1 in
Array.init (Mat.dim1 m_F_mo_a) (fun j ->
let j = j+1 in
match (space i), (space j) with
| `Core , `Core ->
0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) -.
(m_F_mo_b.{i,j} -. m_F_mo_a.{i,j})
| `Active , `Core
| `Core , `Active ->
(*
0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) +.
0.5 *. (m_F_mo_b.{i,j} -. m_F_mo_a.{i,j})
*)
m_F_mo_b.{i,j}
| `Core , `Virtual
| `Virtual , `Core
| `Active , `Active ->
0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j})
| `Virtual , `Active
| `Active , `Virtual ->
(*
0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) -.
0.5 *. (m_F_mo_b.{i,j} -. m_F_mo_a.{i,j})
*)
m_F_mo_a.{i,j}
| `Virtual , `Virtual ->
0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) +.
(m_F_mo_b.{i,j} -. m_F_mo_a.{i,j})
) )
|> Mat.of_array
in
let m_SC =
gemm m_S m_C
in
let m_F =
x_o_xt ~x:m_SC ~o:m_F_mo
in
(* Add level shift in AO basis *)
let m_F =
x_o_xt ~x:m_SC ~o:m_LSmo
|> Mat.add m_F
in
(* 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 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
(* 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_a (Mat.add m_Hc_a m_F_a) +.
Mat.gemm_trace m_P_b (Mat.add m_Hc_b m_F_b) )
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_a m_P_b (Some fock_a) (Some fock_b) error diis
else
let iterations =
List.rev ( (energy, error, gap) :: iterations )
|> Array.of_list
in
HartreeFock_type.(ROHF
{
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_a m_J_a) +.
0.5 *. (Mat.gemm_trace m_P_b m_J_b);
exchange_energy = 0.5 *. (Mat.gemm_trace m_P_a m_K_a) +.
0.5 *. (Mat.gemm_trace m_P_b m_K_b);
occupation = Mat.copy_diag m_P;
})
in
(* Guess coefficients *)
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
let m_C =
gemm m_X m_C'
in
let diis = DIIS.make () in
loop 1 [] None m_C m_C m_C None None threshold_SCF diis

4
Utils/Printing.ml Normal file
View File

@ -0,0 +1,4 @@
let line ?(c='-') n =
String.make n c

4
run_fci.ml
View File

@ -23,9 +23,9 @@ let () =
end;
(* Handle options *)
let basis_file = Util.of_some Command_line.get "basis" in
let basis_file = Util.of_some @@ Command_line.get "basis" in
let nuclei_file = Util.of_some Command_line.get "xyz" in
let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
let charge =
match Command_line.get "charge" with

14
run_hartree_fock.ml
View File

@ -23,9 +23,9 @@ let () =
end;
(* Handle options *)
let basis_file = Util.of_some Command_line.get "basis" in
let basis_file = Util.of_some @@ Command_line.get "basis" in
let nuclei_file = Util.of_some Command_line.get "xyz" in
let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
let charge =
match Command_line.get "charge" with
@ -43,8 +43,12 @@ let () =
Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file
in
HartreeFock.make s
|> HartreeFock.to_string
|> print_endline
let hf =
HartreeFock.make s
in
Format.printf "@[%a@]@." HartreeFock.pp_hf hf;
let mos = MOBasis.of_hartree_fock hf in
Format.printf "@[%a@]@." (fun ppf x -> MOBasis.pp_mo ppf x) mos