mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-08-30 00:03:42 +02:00
F12 almost working
This commit is contained in:
parent
dc50c99b43
commit
29ddf41c0b
2
CI/CI.ml
2
CI/CI.ml
@ -737,7 +737,7 @@ let pt2_en_reference ci =
|
||||
|
||||
let aux_basis = mo_basis in
|
||||
let ds =
|
||||
DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis
|
||||
DeterminantSpace.fci_of_mo_basis ~frozen_core:true aux_basis
|
||||
in
|
||||
let out_dets =
|
||||
ds
|
||||
|
46
CI/F12CI.ml
46
CI/F12CI.ml
@ -136,7 +136,6 @@ Util.debug_matrix "f12 amplitudes" f12_amplitudes;
|
||||
|> Mat.of_array
|
||||
in
|
||||
|
||||
Printf.printf "%d %d\n" (Mat.dim1 m_H_aux) (Mat.dim2 m_H_aux);
|
||||
let m_HF =
|
||||
gemm m_H_aux m_F_aux ~transb:`T
|
||||
in
|
||||
@ -149,7 +148,7 @@ Printf.printf "%d %d\n" (Mat.dim1 m_H_aux) (Mat.dim2 m_H_aux);
|
||||
|
||||
let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =
|
||||
|
||||
let gamma = 1. in
|
||||
let gamma = 1.0 in
|
||||
|
||||
let mo_num = MOBasis.size mo_basis in
|
||||
|
||||
@ -236,9 +235,14 @@ Util.debug_matrix "T" (f12_amplitudes psi);
|
||||
let delta =
|
||||
dressing_vector aux_basis (f12_amplitudes psi) ci
|
||||
in
|
||||
Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
|
||||
(*
|
||||
Util.debug_matrix "psi" psi;
|
||||
*)
|
||||
Format.printf "Amplitude (1,1) : %f@." (f12_amplitudes psi).{1,1};
|
||||
Format.printf "Dressing vector(1,1) : %f@." (Matrix.get delta 1 1);
|
||||
|
||||
(*------
|
||||
TODO SINGLE STATE HERE
|
||||
*)
|
||||
let m_H_dressed = Matrix.to_mat m_H in
|
||||
let f = 1.0 /. psi.{1,1} in
|
||||
@ -251,6 +255,9 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
|
||||
m_H_dressed.{1,i} <- m_H_dressed.{1,i} +. delta;
|
||||
m_H_dressed.{1,1} <- m_H_dressed.{1,1} -. delta *. psi.{i,1} *. f;
|
||||
end
|
||||
(* DIAGONAL DRESSING
|
||||
m_H_dressed.{i,i} <- m_H_dressed.{i,i} +. (Matrix.get delta i 1) /. (psi.{i,1} +. 1.e-5);
|
||||
*)
|
||||
);
|
||||
let eigenvectors, eigenvalues =
|
||||
Util.diagonalize_symm m_H_dressed
|
||||
@ -260,7 +267,7 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
|
||||
(Mat.to_col_vecs psi).(0)
|
||||
(Mat.to_col_vecs eigenvectors).(0) )
|
||||
in
|
||||
Printf.printf "Convergence : %f %f\n" conv (eigenvalues.{1} +. e_shift);
|
||||
Printf.printf "Convergence : %e %f\n" conv (eigenvalues.{1} +. e_shift);
|
||||
if conv > threshold then
|
||||
iteration eigenvectors
|
||||
else
|
||||
@ -272,13 +279,20 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
|
||||
iteration ci_coef
|
||||
(*
|
||||
------- *)
|
||||
|
||||
(*
|
||||
(*TODO : enlever le double comptage de la symmetrisation*)
|
||||
let n_states = ci.CI.n_states in
|
||||
let diagonal =
|
||||
Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i +. (if i=1 then Matrix.get delta 1 1 else 0.) )
|
||||
let f = 1.0 /. (psi.{1,1} ) in
|
||||
let delta_00 =
|
||||
Util.list_range 2 (Mat.dim1 psi)
|
||||
|> List.fold_left (fun accu i -> accu +.
|
||||
(Matrix.get delta i 1) *. psi.{i,1} *. f) 0.
|
||||
in
|
||||
let delta = Matrix.to_mat delta in
|
||||
delta.{1,1} <- delta.{1,1} -. delta_00;
|
||||
let diagonal =
|
||||
Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i +. (if i=1 then delta.{1,1} *. psi.{1,1} else 0.) )
|
||||
in
|
||||
let delta = Matrix.dense_of_mat delta in
|
||||
let matrix_prod c =
|
||||
Matrix.add
|
||||
(Matrix.mm ~transa:`T m_H c)
|
||||
@ -286,24 +300,24 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta;
|
||||
in
|
||||
let eigenvectors, eigenvalues =
|
||||
Parallel.broadcast (lazy (
|
||||
Davidson.make ~threshold:1.e-6 ~guess:(Matrix.to_mat psi) ~n_states diagonal matrix_prod
|
||||
Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
|
||||
))
|
||||
in
|
||||
let m =
|
||||
Matrix.mm ~transa:`T psi (Matrix.dense_of_mat eigenvectors)
|
||||
|> Matrix.to_mat
|
||||
let conv =
|
||||
1.0 -. abs_float ( dot
|
||||
(Mat.to_col_vecs psi).(0)
|
||||
(Mat.to_col_vecs eigenvectors).(0) )
|
||||
in
|
||||
let conv = Mat.sum m -. (Vec.sum (Mat.copy_diag m)) in
|
||||
Printf.printf "Convergence : %f %f\n" conv eigenvalues.{1};
|
||||
Printf.printf "Convergence : %e %f\n" conv (eigenvalues.{1} +. e_shift);
|
||||
if conv > threshold then
|
||||
iteration (Matrix.dense_of_mat eigenvectors)
|
||||
iteration eigenvectors
|
||||
else
|
||||
let eigenvalues =
|
||||
Vec.map (fun x -> x +. e_shift) eigenvalues
|
||||
in
|
||||
eigenvectors, eigenvalues
|
||||
in
|
||||
iteration (Matrix.dense_of_mat ci_coef)
|
||||
iteration ci_coef
|
||||
*)
|
||||
|
||||
)
|
||||
|
@ -56,7 +56,7 @@ let () =
|
||||
in
|
||||
|
||||
let space =
|
||||
DeterminantSpace.fci_of_mo_basis ~frozen_core:false mos
|
||||
DeterminantSpace.fci_of_mo_basis ~frozen_core:false mos
|
||||
in
|
||||
let ci = CI.make space in
|
||||
Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
|
||||
|
Loading…
Reference in New Issue
Block a user