2018-02-23 18:44:31 +01:00
|
|
|
open Util
|
2018-07-04 19:21:45 +02:00
|
|
|
open Constants
|
2018-02-23 18:44:31 +01:00
|
|
|
open Lacaml.D
|
2018-06-13 17:49:58 +02:00
|
|
|
|
|
|
|
module Si = Simulation
|
|
|
|
module El = Electrons
|
|
|
|
module Ao = AOBasis
|
2018-06-13 19:03:42 +02:00
|
|
|
module Ov = Overlap
|
2018-02-23 18:44:31 +01:00
|
|
|
|
2018-05-31 16:46:45 +02:00
|
|
|
let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=0.1)
|
2018-05-31 17:48:54 +02:00
|
|
|
?threshold_SCF:(threshold_SCF=1.e-5) simulation =
|
2018-02-23 18:44:31 +01:00
|
|
|
|
|
|
|
(* Number of occupied MOs *)
|
|
|
|
let nocc =
|
2019-02-20 18:15:15 +01:00
|
|
|
(Si.electrons simulation).El.n_alpha
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
|
|
|
|
2018-05-30 18:07:05 +02:00
|
|
|
let nuclear_repulsion =
|
2019-02-20 18:15:15 +01:00
|
|
|
Si.nuclear_repulsion simulation
|
2018-06-13 17:49:58 +02:00
|
|
|
in
|
|
|
|
|
|
|
|
let ao_basis =
|
2019-02-20 18:15:15 +01:00
|
|
|
Si.ao_basis simulation
|
2018-05-30 18:07:05 +02:00
|
|
|
in
|
|
|
|
|
2018-02-23 18:44:31 +01:00
|
|
|
(* Initial guess *)
|
|
|
|
let guess =
|
2018-06-13 17:49:58 +02:00
|
|
|
Guess.make ~nocc ~guess ao_basis
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
|
|
|
|
|
|
|
(* Orthogonalization matrix *)
|
|
|
|
let m_X =
|
2018-06-13 17:49:58 +02:00
|
|
|
Lazy.force ao_basis.Ao.ortho
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
|
|
|
|
2018-05-28 19:58:40 +02:00
|
|
|
|
2018-02-23 18:44:31 +01:00
|
|
|
(* Overlap matrix *)
|
|
|
|
let m_S =
|
2018-06-13 17:49:58 +02:00
|
|
|
Lazy.force ao_basis.Ao.overlap
|
2018-06-13 19:03:42 +02:00
|
|
|
|> Ov.matrix
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
|
|
|
|
2018-06-13 19:03:42 +02:00
|
|
|
let m_T = Lazy.force ao_basis.Ao.kin_ints |> KinInt.matrix
|
|
|
|
and m_V = Lazy.force ao_basis.Ao.eN_ints |> NucInt.matrix
|
2018-05-30 18:07:05 +02:00
|
|
|
in
|
2018-02-23 18:44:31 +01:00
|
|
|
|
2018-05-31 17:48:54 +02:00
|
|
|
(* Level shift in MO basis *)
|
|
|
|
let m_LSmo =
|
2018-05-31 16:46:45 +02:00
|
|
|
Array.init (Mat.dim2 m_X) (fun i ->
|
|
|
|
if i > nocc then level_shift else 0.)
|
|
|
|
|> Vec.of_array
|
|
|
|
|> Mat.of_diag
|
|
|
|
in
|
|
|
|
|
2018-02-23 18:44:31 +01:00
|
|
|
(* SCF iterations *)
|
2018-07-04 19:21:45 +02:00
|
|
|
let rec loop nSCF iterations energy_prev m_C m_P_prev fock_prev threshold diis =
|
2018-02-23 18:44:31 +01:00
|
|
|
|
|
|
|
(* 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 *)
|
2018-07-04 19:21:45 +02:00
|
|
|
let fock =
|
|
|
|
match fock_prev, threshold > 100. *. threshold_SCF with
|
|
|
|
| Some fock_prev, true ->
|
|
|
|
let threshold = 1.e-8 in
|
|
|
|
Fock.make ~density:(Mat.sub m_P m_P_prev) ~threshold ao_basis
|
|
|
|
|> Fock.add fock_prev
|
|
|
|
| _ -> Fock.make ~density:m_P ao_basis
|
|
|
|
in
|
|
|
|
|
2018-05-30 18:07:05 +02:00
|
|
|
let m_F, m_Hc, m_J, m_K =
|
2018-07-04 19:21:45 +02:00
|
|
|
let x = fock in
|
2018-05-30 18:07:05 +02:00
|
|
|
x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
2018-07-04 19:21:45 +02:00
|
|
|
|
2018-05-31 17:48:54 +02:00
|
|
|
(* 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
|
|
|
|
|
2018-05-31 16:46:45 +02:00
|
|
|
let error_fock =
|
|
|
|
let fps =
|
|
|
|
gemm m_F (gemm m_P m_S)
|
|
|
|
and spf =
|
|
|
|
gemm m_S (gemm m_P m_F)
|
|
|
|
in
|
2018-05-31 17:48:54 +02:00
|
|
|
xt_o_x (Mat.sub fps spf) m_X
|
2018-05-31 16:46:45 +02:00
|
|
|
in
|
|
|
|
|
|
|
|
let diis =
|
2018-05-31 17:48:54 +02:00
|
|
|
DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis
|
2018-05-31 16:46:45 +02:00
|
|
|
in
|
|
|
|
|
|
|
|
let m_F_diis =
|
|
|
|
let x =
|
|
|
|
Bigarray.genarray_of_array1 (DIIS.next diis)
|
|
|
|
in
|
2018-05-31 17:48:54 +02:00
|
|
|
Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho)
|
2018-05-31 16:46:45 +02:00
|
|
|
in
|
|
|
|
|
|
|
|
|
2018-05-31 17:48:54 +02:00
|
|
|
(* MOs in orthogonal MO basis *)
|
2018-02-23 18:44:31 +01:00
|
|
|
let m_C', eigenvalues =
|
2018-05-31 17:48:54 +02:00
|
|
|
diagonalize_symm m_F_diis
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
|
|
|
|
|
|
|
(* MOs in AO basis *)
|
|
|
|
let m_C =
|
2018-05-31 17:48:54 +02:00
|
|
|
gemm m_X m_C'
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
|
|
|
|
|
|
|
(* Hartree-Fock energy *)
|
|
|
|
let energy =
|
2018-05-30 18:07:05 +02:00
|
|
|
nuclear_repulsion +. 0.5 *.
|
2018-02-23 18:44:31 +01:00
|
|
|
Mat.gemm_trace m_P (Mat.add m_Hc m_F)
|
|
|
|
in
|
|
|
|
|
|
|
|
(* Convergence criterion *)
|
2018-05-31 16:46:45 +02:00
|
|
|
let error =
|
2018-05-31 17:48:54 +02:00
|
|
|
error_fock
|
2018-05-31 16:46:45 +02:00
|
|
|
|> Mat.as_vec
|
2018-02-23 18:44:31 +01:00
|
|
|
|> amax
|
2018-05-31 16:46:45 +02:00
|
|
|
|> abs_float
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
|
|
|
|
|
|
|
let converged =
|
2018-05-31 16:46:45 +02:00
|
|
|
nSCF = max_scf || error < threshold_SCF
|
2018-05-30 18:07:05 +02:00
|
|
|
in
|
|
|
|
|
|
|
|
let gap =
|
|
|
|
eigenvalues.{nocc+1} -. eigenvalues.{nocc};
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
|
|
|
|
2018-05-31 17:48:54 +02:00
|
|
|
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
|
2018-02-23 18:44:31 +01:00
|
|
|
|
|
|
|
if not converged then
|
2018-07-04 19:21:45 +02:00
|
|
|
loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C m_P (Some fock) error diis
|
2018-02-23 18:44:31 +01:00
|
|
|
else
|
|
|
|
let iterations =
|
2018-05-31 16:46:45 +02:00
|
|
|
List.rev ( (energy, error, gap) :: iterations )
|
2018-02-23 18:44:31 +01:00
|
|
|
|> Array.of_list
|
|
|
|
in
|
2018-07-20 16:09:06 +02:00
|
|
|
HartreeFock_type.(RHF
|
|
|
|
{
|
|
|
|
simulation;
|
2018-05-30 18:07:05 +02:00
|
|
|
nocc;
|
2018-02-23 18:44:31 +01:00
|
|
|
guess ;
|
|
|
|
eigenvectors = m_C ;
|
|
|
|
eigenvalues ;
|
2018-05-30 18:07:05 +02:00
|
|
|
energy ;
|
|
|
|
nuclear_repulsion;
|
2018-02-23 18:44:31 +01:00
|
|
|
iterations ;
|
2018-05-30 18:07:05 +02:00
|
|
|
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;
|
2018-07-20 16:09:06 +02:00
|
|
|
})
|
2018-02-23 18:44:31 +01:00
|
|
|
in
|
|
|
|
|
|
|
|
|
|
|
|
(* Guess coefficients *)
|
|
|
|
let m_H =
|
|
|
|
match guess with
|
|
|
|
| Guess.Hcore m_H -> m_H
|
2018-05-31 16:46:45 +02:00
|
|
|
| Guess.Huckel m_H -> m_H
|
2018-02-23 18:44:31 +01:00
|
|
|
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
|
|
|
|
|
2018-05-31 16:46:45 +02:00
|
|
|
let diis = DIIS.make () in
|
2018-07-04 19:21:45 +02:00
|
|
|
loop 1 [] None m_C m_C None threshold_SCF diis
|
2018-02-23 18:44:31 +01:00
|
|
|
|
|
|
|
|
|
|
|
|