diff --git a/CI/CI.ml b/CI/CI.ml index d0b57c5..66d2029 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -410,7 +410,7 @@ let make ?(n_states=1) det_space = let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } - list_holes1 list_particles1 + list_holes1 list_particles1 ?(unique=true) list_holes2 list_particles2 i_o1_alfa alfa_o2_i w_alfa psi0 = @@ -462,16 +462,20 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } let det_contribution i = - let already_generated alfa = - if is_internal alfa then - true - else - let rec aux = function - | -1 -> false - | j -> Determinant.degree (fst psi0.(j)) alfa <= 2 - || aux (j-1) - in - aux (i-1) + let already_generated = + if unique then + (fun alfa -> + if is_internal alfa then + true + else + let rec aux = function + | -1 -> false + | j -> Determinant.degree (fst psi0.(j)) alfa <= 2 + || aux (j-1) + in + aux (i-1) + ) else + is_internal in let psi_filtered_idx = diff --git a/CI/F12CI.ml b/CI/F12CI.ml index 4520d20..061835d 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -1,3 +1,6 @@ +let debug s = + Printf.printf "%s\n%!" s; + open Lacaml.D type t = @@ -7,13 +10,13 @@ type t = det_space : DeterminantSpace.t ; ci : CI.t ; eigensystem : (Mat.t * Vec.t) lazy_t; - f12_amplitudes : Mat.t; } let ci t = t.ci let mo_basis t = t.mo_basis let det_space t = t.det_space let mo_class t = DeterminantSpace.mo_class @@ det_space t +let eigensystem t = Lazy.force t.eigensystem let f12_integrals mo_basis = @@ -50,6 +53,7 @@ let f_ij mo_basis ki kj = let dressing_vector f12_amplitudes ci = +debug "Computing dressing vector"; let mo_basis = DeterminantSpace.mo_basis ci.CI.det_space in let i_o1_alfa = h_ij mo_basis in @@ -65,10 +69,11 @@ let dressing_vector f12_amplitudes ci = and list_particles2 = List.concat [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ; MOClass.auxiliary_mos mo_class ] in +Util.debug_matrix "f12 amplitudes" f12_amplitudes; (* 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 + i_o1_alfa alfa_o2_i w_alfa f12_amplitudes ~unique:false |> Vec.of_list in Matrix.sparse_of_vector_array [| Vector.sparse_of_vec result |] @@ -108,51 +113,94 @@ let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basi let ci = CI.make det_space in - let ci_coef, ci_energy = Parallel.broadcast ci.eigensystem in + let ci_coef, ci_energy = + let x = Lazy.force ci.eigensystem in + Parallel.broadcast (lazy x) + in let f12_amplitudes = (* While in a sequential region, initiate the parallel 4-idx transformation to avoid nested parallel jobs *) - ignore @@ MOBasis.f12_ints mo_basis; +debug "Four-idx transform of f12 intergals"; + ignore @@ MOBasis.f12_ints aux_basis; let f = fun ki kj -> if ki <> kj then - f_ij mo_basis ki kj + f_ij aux_basis ki kj else - f_ij mo_basis ki kj +. 1. + f_ij aux_basis ki kj +. 1. in +debug "Computing F matrix"; let m_F = CI.create_matrix_spin f det_space |> Lazy.force in + fun ci_coef -> +debug "Solving linear system"; Matrix.ax_eq_b m_F (Matrix.dense_of_mat ci_coef) |> Matrix.to_mat in - let f_11, e_shift = + let e_shift = let det = DeterminantSpace.determinant_stream det_space |> Stream.next in - f_ij mo_basis det det, - h_ij mo_basis det det + h_ij aux_basis det det in 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. ) +debug "Iteration"; + let delta = + dressing_vector (f12_amplitudes psi) ci in - let matrix_prod psi = +Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta; + +(*------ +(*TODO : enlever le double comptage de la symmetrisation*) +*) + let m_H_dressed = Matrix.to_mat m_H in + (Matrix.dim1 delta) (Matrix.dim2 delta) (Mat.dim1 psi) (Mat.dim2 psi) ; + Util.list_range 1 (Mat.dim1 psi) + |> List.iter (fun i -> m_H_dressed.{i,1} <- m_H_dressed.{i,1} +. (Matrix.get delta i 1) /. (psi.{1,1})); + let eigenvectors, eigenvalues = + Util.diagonalize_symm m_H_dressed + in + let m = + gemm ~transa:`T psi eigenvectors + |> Mat.abs + 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 eigenvectors + else + let eigenvalues = + Vec.map (fun x -> x +. e_shift) eigenvalues + in + eigenvectors, eigenvalues + in + iteration ci_coef +(* +------- *) + +(* + 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.) ) + in + let matrix_prod c = Matrix.add - (Matrix.mm ~transa:`T m_H psi) - (dressing_vector f12_amplitudes ci) + (Matrix.mm ~transa:`T m_H c) + delta in let eigenvectors, eigenvalues = Parallel.broadcast (lazy ( @@ -174,8 +222,10 @@ let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basi eigenvectors, eigenvalues in iteration (Matrix.dense_of_mat ci_coef) +*) + ) in - { mo_basis ; aux_basis ; det_space ; ci ; f12_amplitudes ; eigensystem } + { mo_basis ; aux_basis ; det_space ; ci ; eigensystem } diff --git a/Utils/Matrix.ml b/Utils/Matrix.ml index 96c0164..9ff0514 100644 --- a/Utils/Matrix.ml +++ b/Utils/Matrix.ml @@ -35,7 +35,7 @@ let dim2 = function let get = function | Dense m -> (fun i j -> m.{i,j}) - | Sparse {m ; n ; v } -> (fun i j -> Vector.get v.(i-1) j) + | Sparse {m ; n ; v } -> (fun i j -> Vector.get v.(j-1) i) let sparse_of_dense ?(threshold=epsilon) = function diff --git a/run_fci_f12.ml b/run_fci_f12.ml index ddc20f1..08e5325 100644 --- a/run_fci_f12.ml +++ b/run_fci_f12.ml @@ -65,10 +65,6 @@ let () = F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename () in let ci = F12CI.ci fcif12 in - Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion simulation); - (* - let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in - Util.list_range 1 (DeterminantSpace.size space) - |> List.iter (fun i -> Format.printf "@[%f@]@;" s2.{i,i}); - *) - + Format.fprintf ppf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion simulation); + let _, e_cif12 = F12CI.eigensystem fcif12 in + Format.fprintf ppf "FCI-F12 energy : %20.16f@." (e_cif12.{1} +. Simulation.nuclear_repulsion simulation);