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

Working on f12.

This commit is contained in:
Anthony Scemama 2019-05-07 17:00:44 +02:00
parent dea45d7b52
commit d1efe67e5a
5 changed files with 98 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Basis/F12.ml
View File

@ -373,9 +373,6 @@ let of_basis_parallel f12 basis =
Fis.create ~size:n `Dense
in
(*
let lambda_inv = -. 1. /. f12.expo_s in
*)
Farm.run ~ordered:false ~f input_stream
|> Stream.iter (fun l ->
Array.iter (fun (i_c,j_c,k_c,l_c,value) ->

53
CI/F12CI.ml
View File

@ -24,11 +24,16 @@ let f12_integrals mo_basis =
0.
else
begin
let ijkl = F12.get_phys two_e_ints i j k l
and ijlk = F12.get_phys two_e_ints i j l k
in
if s' = Spin.other s then
0.5 *. (F12.get_phys two_e_ints i j k l)
0.375 *. ijkl +. 0.125 *. ijlk
(*
0.25 *. ijkl
*)
else
0.25 *. ((F12.get_phys two_e_ints i j k l) -.
(F12.get_phys two_e_ints i j l k))
0.25 *. (ijkl -. ijlk)
end
) )
@ -204,7 +209,7 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci =
let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () =
let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename ?(state=1) () =
let f12 = Util.of_some @@ Simulation.f12 simulation in
let mo_num = MOBasis.size mo_basis in
@ -243,7 +248,7 @@ Printf.printf "det space\n%!";
DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num
in
let ci = CI.make det_space in
let ci = CI.make ~n_states:state det_space in
let ci_coef, ci_energy =
let x = Lazy.force ci.eigensystem in
@ -263,28 +268,26 @@ Printf.printf "det space\n%!";
Lazy.force ci.CI.m_H
in
let rec iteration ?(state=1) psi =
let rec iteration ~state psi =
let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in
let delta =
dressing_vector ~frozen_core aux_basis psi ci
in
let f = 1.0 /. psi.{1,1} in
let f = 1.0 /. psi.{column_idx,state} 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.
(Matrix.get delta i state) *. psi.{i,state} *. f) 0.
in
let delta = Matrix.to_mat delta in
delta.{1,1} <- delta.{1,1} -. delta_00;
delta.{column_idx,state} <- delta.{column_idx,state} -. delta_00;
(*------
TODO SINGLE STATE HERE
*)
let n_states = ci.CI.n_states in
let diagonal =
Vec.init (Matrix.dim1 m_H) (fun i ->
if i = 1 then
Matrix.get m_H i i +. delta.{1,1} *. f
if i = column_idx then
Matrix.get m_H i i +. delta.{column_idx,state} *. f
else
Matrix.get m_H i i
)
@ -295,24 +298,26 @@ TODO SINGLE STATE HERE
Matrix.mm ~transa:`T m_H c
|> Matrix.to_mat
in
let c11 = Matrix.get c 1 1 in
let c11 = Matrix.get c state state in
Util.list_range 1 (Mat.dim1 w)
|> List.iter (fun i ->
let dci =
delta.{i,1} *. f ;
delta.{i,state} *. f ;
in
w.{i,1} <- w.{i,1} +. dci *. c11;
if (i <> 1) then
w.{1,1} <- w.{1,1} +. dci *. (Matrix.get c i 1);
w.{i,state} <- w.{i,state} +. dci *. c11;
if (i <> column_idx) then
w.{column_idx,state} <- w.{column_idx,state} +. dci *. (Matrix.get c i state);
);
Matrix.dense_of_mat w
in
let eigenvectors, eigenvalues =
Parallel.broadcast (lazy (
Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states diagonal matrix_prod
Davidson.make ~threshold:1.e-6 ~guess:psi ~n_states:state diagonal matrix_prod
))
in
Vec.iter (fun energy -> Printf.printf "%g\t" energy) eigenvalues;
print_newline ();
let conv =
1.0 -. abs_float ( dot
@ -320,18 +325,18 @@ TODO SINGLE STATE HERE
(Mat.to_col_vecs eigenvectors).(0) )
in
if Parallel.master then
Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{1} +. e_shift
Printf.printf "F12 Convergence : %e %f\n" conv (eigenvalues.{state} +. e_shift
+. Simulation.nuclear_repulsion simulation);
if conv > threshold then
iteration eigenvectors
iteration ~state eigenvectors
else
let eigenvalues =
Vec.map (fun x -> x +. e_shift) eigenvalues
in
eigenvectors, eigenvalues
in
iteration ci_coef
iteration ~state ci_coef
)
in

66
Utils/Davidson.ml
View File

@ -14,15 +14,10 @@ let make
(* Size of the matrix to diagonalize *)
let n = Vec.dim diagonal in
let m = (* Number of requested states *)
match guess with
| Some vectors -> (Mat.dim2 vectors) * n_states
| None -> n_states
in
let m = n_states in
(* Create guess vectors u, with randomly initialized unknown vectors. *)
let init_vectors =
let init_vectors m =
let init_vector k =
Vector.sparse_of_assoc_list n [ (k,1.0) ]
in
@ -31,6 +26,20 @@ let make
|> Matrix.to_mat
in
let guess =
match guess with
| Some vectors -> vectors
| None -> init_vectors m
in
let guess =
if Mat.dim2 guess = n_states then guess
else
(Mat.to_col_vecs_list guess) @
(Mat.to_col_vecs_list (init_vectors (m-(Mat.dim2 guess))) )
|> Mat.of_col_vecs_list
in
let pick_new u =
Mat.to_col_vecs_list u
|> Util.list_pack m
@ -38,19 +47,14 @@ let make
|> List.hd
in
let u_new =
match guess with
| Some vectors -> Mat.to_col_vecs_list vectors
| None -> Mat.to_col_vecs_list init_vectors
in
let u_new = Mat.to_col_vecs_list guess in
let rec iteration u u_new w iter =
let rec iteration u u_new w iter macro =
(* u is a list of orthonormal vectors, on which the operator has
been applied : w = op.u
u_new is a list of vector which will increase the size of the
u_new is a list of vectors which will increase the size of the
space.
*)
(* Orthonormalize input vectors u_new *)
let u_new_ortho =
List.concat [u ; u_new]
@ -101,17 +105,29 @@ let make
(* Compute the residual as proposed new vectors *)
let u_proposed =
Mat.init_cols n m (fun i k ->
(lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /.
(max (diagonal.{i} -. lambda.{k}) 0.01) )
let delta = diagonal.{i} -. lambda.{k} in
let delta =
if abs_float delta > 0.01 then delta
else if delta < 0. then -.0.01
else 0.01
in
(lambda.{k} *. m_new_U.{i,k} -. m_new_W.{i,k}) /. delta )
|> Mat.to_col_vecs_list
in
let residual_norms = List.map nrm2 u_proposed in
let residual_norm = List.fold_left (fun accu i -> max accu i) 0. residual_norms in
let residual_norm =
List.fold_left (fun accu i -> accu +. i *. i) 0. residual_norms
|> sqrt
in
if Parallel.master then
Printf.printf "%3d %16.10f %16.8e%!\n" iter lambda.{1} residual_norm;
begin
Printf.printf "%3d " iter;
Vec.iteri (fun i x -> if (i<=m) then Printf.printf "%16.10f " x) lambda;
Printf.printf "%16.8e%!\n" residual_norm
end;
(* Make new vectors sparse *)
let u_proposed =
@ -125,20 +141,20 @@ let make
in
if residual_norm > threshold then
let u_next, w_next, iter =
if residual_norm > threshold && macro > 0 then
let u_next, w_next, iter, macro =
if iter = n_iter then
m_new_U |> pick_new,
m_new_W |> pick_new,
0
0, (macro-1)
else
u_next, w_next, iter
u_next, w_next, (iter+1), macro
in
iteration u_next u_proposed w_next (iter+1)
iteration u_next u_proposed w_next iter macro
else
(m_new_U |> pick_new |> Mat.of_col_vecs_list), lambda
in
iteration [] u_new [] 1
iteration [] u_new [] 1 5

28
run_fci_f12.ml
View File

@ -1,3 +1,5 @@
open Lacaml.D
let () =
let open Command_line in
begin
@ -27,6 +29,10 @@ let () =
{ short='e' ; long="expo" ; opt=Optional;
arg=With_arg "<float>";
doc="Exponent of the Gaussian geminal"; } ;
{ short='s' ; long="state" ; opt=Optional;
arg=With_arg "<int>";
doc="Requested state. Default is 1."; } ;
]
end;
@ -49,6 +55,12 @@ let () =
| None -> 1.0
in
let state =
match Command_line.get "state" with
| Some x -> int_of_string x
| None -> 1
in
let multiplicity =
match Command_line.get "multiplicity" with
| Some x -> int_of_string x
@ -76,11 +88,21 @@ let () =
in
let fcif12 =
F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ()
F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ~state ()
in
let ci = F12CI.ci fcif12 in
Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion simulation);
Format.fprintf ppf "FCI energy : ";
Vec.iteri (fun i x -> if i <= state then
Format.fprintf ppf "%20.16f@; " (x +. Simulation.nuclear_repulsion simulation) )
(CI.eigenvalues ci);
Format.fprintf ppf "@.";
let _, e_cif12 = F12CI.eigensystem fcif12 in
Format.fprintf ppf "FCI-F12 energy : %20.16f@." (e_cif12.{1} +. Simulation.nuclear_repulsion simulation);
Format.fprintf ppf "FCI-F12 energy : ";
Vec.iteri (fun i x -> if i <= state then
Format.fprintf ppf "%20.16f@; " (x +. Simulation.nuclear_repulsion simulation) )
e_cif12;
Format.fprintf ppf "@."

7
run_integrals.ml
View File

@ -27,8 +27,9 @@ let run ~out =
| Some x -> x
in
let f12 = F12.gaussian_geminal 1.0 in
let s =
Simulation.of_filenames ~nuclei:nuclei_file basis_file
Simulation.of_filenames ~nuclei:nuclei_file ~f12 basis_file
in
print_endline @@ Nuclei.to_string @@ Simulation.nuclei s;
@ -45,10 +46,8 @@ let run ~out =
NucInt.to_file ~filename:(out_file^".nuc") eN_ints;
KinInt.to_file ~filename:(out_file^".kin") kin_ints;
ERI.to_file ~filename:(out_file^".eri") ee_ints;
(*
let f12_ints = AOBasis.f12_ints ao_basis in
let f12_ints = AOBasis.f12_ints ao_basis in
F12.to_file ~filename:(out_file^".f12") f12_ints;
*)
()