diff --git a/Basis/TwoElectronIntegralsSeparable.ml b/Basis/TwoElectronIntegralsSeparable.ml index ebc54cb..b5fe2b8 100644 --- a/Basis/TwoElectronIntegralsSeparable.ml +++ b/Basis/TwoElectronIntegralsSeparable.ml @@ -17,7 +17,7 @@ include FourIdxStorage (** Exponent of the geminal *) let expo_s = 1.0 -(** Coefficients and exponents of the Gaussian fit *) +(** Coefficients and exponents of the Gaussian fit of the Slater Geminal*) let coef_g = [| 0.3144 ; 0.3037 ; 0.1681 ; 0.09811 ; 0.06024 ; 0.03726 |] @@ -26,6 +26,29 @@ let expo_sg_inv = [| 0.2209 ; 1.004 ; 3.622 ; 12.16 ; 45.87 ; 254.4 |] + +(* + +Fit of 1/r: + +let coef_g = [| + 841.88478132 ; 70.590185207 ; 18.3616020768 ; 7.2608642093 ; + 3.57483416444 ; 2.01376031082 ; 1.24216542801 ; 0.81754348620 ; + 0.564546514023 ; 0.404228610699 ; 0.297458536575 ; 0.223321219537 ; + 0.169933732064 ; 0.130190978230 ; 0.099652303426 ; 0.075428246546 ; + 0.0555635614051 ; 0.0386791283055 ; 0.0237550435652 ; 0.010006278387 ; + |] + +let expo_sg_inv = + Array.map (fun x -> 1. /. (x *. expo_s *. expo_s)) + [| 84135.654509 ; 2971.58727634 ; 474.716025959 ; 130.676724560 ; + 47.3938388887 ; 20.2078651631 ; 9.5411021938 ; 4.8109546955 ; + 2.52795733067 ; 1.35894103210 ; 0.73586710268 ; 0.39557629706 ; + 0.20785895177 ; 0.104809693858 ; 0.049485682527 ; 0.021099788990 ; + 0.007652472186 ; 0.0021065225215 ; 0.0003365204879 ; 0.0000118855674 |] +*) + + let class_of_contracted_shell_pair_couple shell_pair_couple = F12RR.contracted_class_shell_pair_couple expo_sg_inv coef_g shell_pair_couple diff --git a/CI/CI.ml b/CI/CI.ml index ac7d336..88ffd91 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -5,6 +5,7 @@ module Sd = Spindeterminant type t = { + e_shift : float ; (* Diagonal energy shift for increasing numerical precision *) det_space : Ds.t ; m_H : Matrix.t lazy_t ; m_S2 : Matrix.t lazy_t ; @@ -342,9 +343,17 @@ let create_matrix_spin f det_space = let make ?(n_states=1) det_space = - let m_H = + let mo_basis = Ds.mo_basis det_space in - let mo_basis = Ds.mo_basis det_space in + let e_shift = + let d0 = + Ds.determinant_stream det_space + |> Stream.next + in + h_ij mo_basis d0 d0 + in + + let m_H = (* While in a sequential region, initiate the parallel 4-idx transformation to avoid nested parallel jobs @@ -356,7 +365,12 @@ let make ?(n_states=1) det_space = | Ds.Arbitrary _ -> create_matrix_arbitrary | Ds.Spin _ -> create_matrix_spin in - f (fun ki kj -> h_ij mo_basis ki kj) det_space + f (fun ki kj -> + if ki <> kj then + h_ij mo_basis ki kj + else + h_ij mo_basis ki kj -. e_shift + ) det_space in let m_S2 = @@ -378,10 +392,14 @@ let make ?(n_states=1) det_space = let matrix_prod psi = Matrix.mm ~transa:`T m_H psi in - Davidson.make ~n_states diagonal matrix_prod + let eigenvectors, eigenvalues = + Davidson.make ~threshold:1.e-12 ~n_states diagonal matrix_prod + in + let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in + eigenvectors, eigenvalues ) in - { det_space ; m_H ; m_S2 ; eigensystem ; n_states } + { det_space ; e_shift ; m_H ; m_S2 ; eigensystem ; n_states } @@ -411,8 +429,6 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } Array.init (Ds.size det_space) (fun i -> Stream.next stream, psi0.{i+1,1}) in - - (* let is_internal = let m l = @@ -447,8 +463,8 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let is_internal alfa = let rec aux = function | -1 -> false - | j -> if (Determinant.degree (fst psi0.(j)) alfa = 0) then true - else aux (j-1) + | j -> Determinant.degree (fst psi0.(j)) alfa = 0 + || aux (j-1) in aux (Array.length psi0 - 1) in @@ -462,8 +478,8 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } else let rec aux = function | -1 -> false - | j -> if (Determinant.degree (fst psi0.(j)) alfa <= 2) then true - else aux (j-1) + | j -> Determinant.degree (fst psi0.(j)) alfa <= 2 + || aux (j-1) in aux (i-1) in @@ -506,6 +522,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let det_i = fst psi0.(i) in let w_alfa = w_alfa det_i in + let same_spin = List.fold_left (fun accu spin -> accu +. @@ -525,12 +542,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let double = List.fold_left (fun accu particle' -> - if particle' > particle then + if particle' > particle || particle' = hole then accu else accu +. List.fold_left (fun accu hole' -> - if hole' = particle' || hole' < hole then + if hole' = particle' || hole' = particle || hole' < hole then accu else let alfa = diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml index d0e078f..a6b313a 100644 --- a/Utils/Davidson.ml +++ b/Utils/Davidson.ml @@ -11,7 +11,8 @@ let make matrix_prod = - let n = Vec.dim diagonal in (* Size of the matrix to diagonalize *) + (* Size of the matrix to diagonalize *) + let n = Vec.dim diagonal in let m = (* Number of requested states *) @@ -27,6 +28,8 @@ let make if i 1.e-1 then c else 0. + *) ) |> Util.normalize in @@ -130,7 +134,7 @@ let make in iteration u_next u_proposed w_next (iter+1) else - (Mat.of_col_vecs_list u_next |> pick_new |> Mat.of_col_vecs_list), lambda + (m_new_U |> pick_new |> Mat.of_col_vecs_list), lambda in iteration [] u_new [] 1