From d1efe67e5ae2968ad8bd6e6c75b184413937f3f5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 7 May 2019 17:00:44 +0200 Subject: [PATCH] Working on f12. --- Basis/F12.ml | 3 --- CI/F12CI.ml | 53 ++++++++++++++++++++----------------- Utils/Davidson.ml | 66 +++++++++++++++++++++++++++++------------------ run_fci_f12.ml | 28 +++++++++++++++++--- run_integrals.ml | 7 +++-- 5 files changed, 98 insertions(+), 59 deletions(-) diff --git a/Basis/F12.ml b/Basis/F12.ml index 2540ea8..3f0c79f 100644 --- a/Basis/F12.ml +++ b/Basis/F12.ml @@ -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) -> diff --git a/CI/F12CI.ml b/CI/F12CI.ml index f76f8d2..9beb0f4 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -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 diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml index bd52d7f..d72d482 100644 --- a/Utils/Davidson.ml +++ b/Utils/Davidson.ml @@ -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 diff --git a/run_fci_f12.ml b/run_fci_f12.ml index 5a374a0..d7e3390 100644 --- a/run_fci_f12.ml +++ b/run_fci_f12.ml @@ -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 ""; doc="Exponent of the Gaussian geminal"; } ; + + { short='s' ; long="state" ; opt=Optional; + arg=With_arg ""; + 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 "@." + + diff --git a/run_integrals.ml b/run_integrals.ml index 5b69c30..98dd70e 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -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; - *) ()