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

Huckel guess and DIIS OK

This commit is contained in:
Anthony Scemama 2018-05-31 17:48:54 +02:00
parent 3482e02695
commit e81a48dc73
2 changed files with 36 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
SCF/RHF.ml
View File

@ -3,7 +3,7 @@ open Lacaml.D
open Simulation
let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=0.1)
?threshold_SCF:(threshold_SCF=1.e-6) simulation =
?threshold_SCF:(threshold_SCF=1.e-5) simulation =
(* Number of occupied MOs *)
let nocc =
@ -34,8 +34,8 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
and m_V = Lazy.force simulation.eN_ints
in
(* Level shift *)
let m_LS =
(* 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
@ -43,7 +43,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
in
(* SCF iterations *)
let rec loop nSCF iterations m_C diis =
let rec loop nSCF iterations energy_prev m_C diis =
(* Density matrix over nocc occupied MOs *)
let m_P =
@ -58,41 +58,50 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange
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
Mat.sub fps spf
xt_o_x (Mat.sub fps spf) m_X
in
let diis =
DIIS.append ~p:(Mat.as_vec m_F) ~e:(Mat.as_vec error_fock) 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) (Mat.dim2 m_F)
Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho)
in
(* Fock matrix in MO basis *)
let m_Fmo =
xt_o_x m_F_diis m_C
|> Mat.add m_LS
in
(* MOs in old MO basis *)
(* MOs in orthogonal MO basis *)
let m_C', eigenvalues =
diagonalize_symm m_Fmo
diagonalize_symm m_F_diis
in
(* MOs in AO basis *)
let m_C =
gemm m_C m_C'
gemm m_X m_C'
in
(* Hartree-Fock energy *)
@ -103,7 +112,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
(* Convergence criterion *)
let error =
xt_o_x error_fock m_C
error_fock
|> Mat.as_vec
|> amax
|> abs_float
@ -117,10 +126,16 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
eigenvalues.{nocc+1} -. eigenvalues.{nocc};
in
Printf.printf "%d %16.10f %11.4e %10.4f\n%!" nSCF energy error gap;
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) m_C diis
loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C diis
else
let iterations =
List.rev ( (energy, error, gap) :: iterations )
@ -159,7 +174,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=
in
let diis = DIIS.make () in
loop 1 [] m_C diis
loop 1 [] None m_C diis

2
Utils/DIIS.ml
View File

@ -44,7 +44,7 @@ let next diis =
let a = Mat.make (m+1) (m+1) 1. in
a.{m+1,m+1} <- 0.;
ignore @@ lacpy ~b:a (gemm ~transa:`T ~m ~n:m e e);
if sycon (lacpy a) > 1.e-10 then a
if sycon (lacpy a) > 1.e-14 then a
else aux (m-1)
in
aux diis.m