From 0d935ea35484a8d4b63b6b465dd68394c7fcdbdb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 21 Mar 2019 21:48:21 +0100 Subject: [PATCH] Not working yet --- CI/CI.ml | 40 +++++++++++++---------- CI/F12CI.ml | 85 +++++++++++++++++++++++++++++++++++++++--------- Utils/Matrix.mli | 5 +++ 3 files changed, 97 insertions(+), 33 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index f8f6ec4..d0b57c5 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -410,10 +410,14 @@ let make ?(n_states=1) det_space = let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } - list_holes list_particles i_o1_alfa alfa_o2_i w_alfa psi0 = + list_holes1 list_particles1 + list_holes2 list_particles2 + i_o1_alfa alfa_o2_i w_alfa psi0 = - let list_holes = Array.of_list list_holes - and list_particles = Array.of_list list_particles + let list_holes1 = Array.of_list list_holes1 + and list_holes2 = Array.of_list list_holes2 + and list_particles2 = Array.of_list list_particles1 + and list_particles1 = Array.of_list list_particles2 in let psi0 = @@ -422,7 +426,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } Ds.determinant_stream det_space in Array.init (Ds.size det_space) (fun i -> - Stream.next stream, psi0.{i+1,1}) + (Stream.next stream), (Mat.copy_row psi0 (i+1)) ) in let is_internal = @@ -488,7 +492,8 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let psi_h_alfa alfa = List.fold_left (fun accu (det, coef) -> - accu +. coef *. (i_o1_alfa det alfa)) 0. psi_filtered +(* Single state here *) + accu +. coef.{1} *. (i_o1_alfa det alfa)) 0. psi_filtered in let alfa_h_psi = @@ -497,7 +502,8 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } else fun alfa -> List.fold_left (fun accu (det, coef) -> - accu +. coef *. (alfa_o2_i alfa det)) 0. psi_filtered +(* Single state here *) + accu +. coef.{1} *. (alfa_o2_i alfa det)) 0. psi_filtered in let psi_h_alfa_alfa_h_psi alfa = @@ -547,12 +553,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } accu else accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa - ) 0. list_holes - ) 0. list_particles + ) 0. list_holes1 + ) 0. list_particles1 in accu +. single +. double - ) 0. list_holes - ) 0. list_particles + ) 0. list_holes2 + ) 0. list_particles2 ) 0. [ Spin.Alfa ; Spin.Beta ] in @@ -581,12 +587,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } accu else accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa - ) 0. list_holes - ) 0. list_particles + ) 0. list_holes1 + ) 0. list_particles1 in accu +. double - ) 0. list_holes - ) 0. list_particles + ) 0. list_holes2 + ) 0. list_particles2 in same_spin +. opposite_spin in @@ -656,7 +662,7 @@ let pt2_en ci = [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ] in - second_order_sum ci list_holes list_particles + second_order_sum ci list_holes list_particles list_holes list_particles i_o1_alfa i_o1_alfa w_alfa psi0 |> List.fold_left (+.) 0. @@ -688,7 +694,7 @@ let pt2_mp ci = in let psi0, _ = Parallel.broadcast ci.eigensystem in - second_order_sum ci list_holes list_particles + second_order_sum ci list_holes list_particles list_holes list_particles i_o1_alfa i_o1_alfa w_alfa psi0 |> List.fold_left (+.) 0. @@ -709,7 +715,7 @@ let variance ci = in let psi0, _ = Parallel.broadcast ci.eigensystem in - second_order_sum ci list_holes list_particles + second_order_sum ci list_holes list_particles list_holes list_particles i_o1_alfa i_o1_alfa w_alfa psi0 |> List.fold_left (+.) 0. diff --git a/CI/F12CI.ml b/CI/F12CI.ml index e27cd82..4520d20 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -6,6 +6,7 @@ type t = aux_basis : MOBasis.t ; det_space : DeterminantSpace.t ; ci : CI.t ; + eigensystem : (Mat.t * Vec.t) lazy_t; f12_amplitudes : Mat.t; } @@ -47,7 +48,7 @@ let f_ij mo_basis ki kj = |> List.hd -let dressing_vector ci f12_amplitudes = +let dressing_vector f12_amplitudes ci = let mo_basis = DeterminantSpace.mo_basis ci.CI.det_space in @@ -60,17 +61,24 @@ let dressing_vector ci f12_amplitudes = let list_holes = List.concat [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ] - and list_particles = MOClass.auxiliary_mos mo_class + and list_particles1 = MOClass.auxiliary_mos mo_class + and list_particles2 = List.concat + [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ; MOClass.auxiliary_mos mo_class ] in - CI.second_order_sum ci list_holes list_particles - i_o1_alfa alfa_o2_i w_alfa f12_amplitudes - |> Vec.of_list + (* Single state here *) + let result = + CI.second_order_sum ci list_holes list_particles1 list_holes list_particles2 + i_o1_alfa alfa_o2_i w_alfa f12_amplitudes + |> Vec.of_list + in + Matrix.sparse_of_vector_array [| Vector.sparse_of_vec result |] -let make ~simulation ?(frozen_core=true) ~mo_basis ~aux_basis_filename () = + +let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basis_filename () = let mo_num = MOBasis.size mo_basis in @@ -81,11 +89,11 @@ let make ~simulation ?(frozen_core=true) ~mo_basis ~aux_basis_filename () = and nuclei = Simulation.nuclei simulation in let general_basis = - Basis.general_basis @@ Simulation.basis simulation + Basis.general_basis @@ Simulation.basis simulation in GeneralBasis.combine [ - general_basis ; GeneralBasis.read aux_basis_filename - ] + general_basis ; GeneralBasis.read aux_basis_filename + ] |> Basis.of_nuclei_and_general_basis nuclei |> Simulation.make ~charge ~multiplicity ~nuclei in @@ -95,7 +103,7 @@ let make ~simulation ?(frozen_core=true) ~mo_basis ~aux_basis_filename () = in let det_space = - DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num + DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num in let ci = CI.make det_space in @@ -109,20 +117,65 @@ let make ~simulation ?(frozen_core=true) ~mo_basis ~aux_basis_filename () = ignore @@ MOBasis.f12_ints mo_basis; let f = fun ki kj -> - if ki <> kj then - f_ij mo_basis ki kj - else - f_ij mo_basis ki kj +. 1. + if ki <> kj then + f_ij mo_basis ki kj + else + f_ij mo_basis ki kj +. 1. in let m_F = CI.create_matrix_spin f det_space |> Lazy.force in - Matrix.ax_eq_b (Matrix.dense_of_sparse m_F) (Matrix.dense_of_mat ci_coef) + Matrix.ax_eq_b m_F (Matrix.dense_of_mat ci_coef) |> Matrix.to_mat in + let f_11, e_shift = + let det = + DeterminantSpace.determinant_stream det_space + |> Stream.next + in + f_ij mo_basis det det, + h_ij mo_basis det det + in - { mo_basis ; aux_basis ; det_space ; ci ; f12_amplitudes } + let eigensystem = lazy ( + let m_H = + Lazy.force ci.CI.m_H + in + let n_states = ci.CI.n_states in + + let rec iteration ?(state=1) psi = + let diagonal = + Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i +. if i=1 then f_11 /. (Matrix.get psi 1 state) else 0. ) + in + let matrix_prod psi = + Matrix.add + (Matrix.mm ~transa:`T m_H psi) + (dressing_vector f12_amplitudes ci) + in + let eigenvectors, eigenvalues = + Parallel.broadcast (lazy ( + Davidson.make ~threshold:1.e-6 ~guess:(Matrix.to_mat psi) ~n_states diagonal matrix_prod + )) + in + let m = + Matrix.mm ~transa:`T psi (Matrix.dense_of_mat eigenvectors) + |> Matrix.to_mat + in + let conv = Mat.sum m -. (Vec.sum (Mat.copy_diag m)) in + Printf.printf "Convergence : %f %f\n" conv eigenvalues.{1}; + if conv > threshold then + iteration (Matrix.dense_of_mat eigenvectors) + else + let eigenvalues = + Vec.map (fun x -> x +. e_shift) eigenvalues + in + eigenvectors, eigenvalues + in + iteration (Matrix.dense_of_mat ci_coef) + ) + in + { mo_basis ; aux_basis ; det_space ; ci ; f12_amplitudes ; eigensystem } diff --git a/Utils/Matrix.mli b/Utils/Matrix.mli index e743a9f..1bffbc3 100644 --- a/Utils/Matrix.mli +++ b/Utils/Matrix.mli @@ -63,6 +63,11 @@ val mv : ?sparse:bool -> ?trans:trans3 -> ?threshold:float -> t -> Vector.t -> V val ax_eq_b : ?trans:trans3 -> t -> t -> t (** Solves A.X = B or A'.X = B *) +val add : t -> t -> t +(** Add two matrices *) + +val sub : t -> t -> t +(** Subtract two matrices *) (** {1 Printers } *)