From dc50c99b434b8d9c1d56adf58106007def5077fc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 22 Mar 2019 18:16:41 +0100 Subject: [PATCH] Introduced reference PT2 energy --- CI/CI.ml | 124 +++++++++++++++++++++++++++----------- CI/F12CI.ml | 130 ++++++++++++++++++++++++++++++++-------- Utils/FourIdxStorage.ml | 5 ++ run_cas.ml | 10 +++- 4 files changed, 210 insertions(+), 59 deletions(-) diff --git a/CI/CI.ml b/CI/CI.ml index 66d2029..1476f47 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -411,13 +411,13 @@ let make ?(n_states=1) det_space = let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } list_holes1 list_particles1 ?(unique=true) - list_holes2 list_particles2 + list_holes2 list_particles2 is_internal i_o1_alfa alfa_o2_i w_alfa psi0 = 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 + and list_particles1 = Array.of_list list_particles1 + and list_particles2 = Array.of_list list_particles2 in let psi0 = @@ -429,32 +429,6 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } (Stream.next stream), (Mat.copy_row psi0 (i+1)) ) in - let is_internal = - let m l = - List.fold_left (fun accu i -> - let j = i-1 in Z.(logor accu (shift_left one j)) - ) Z.zero l - in - let mo_class = DeterminantSpace.mo_class det_space in - let active_mask = m (MOClass.active_mos mo_class) in - let occ_mask = m (MOClass.core_mos mo_class) in - let inactive_mask = m (MOClass.inactive_mos mo_class) in - let occ_mask = Z.logor occ_mask inactive_mask in - let neg_active_mask = Z.lognot active_mask in - fun a -> - let alfa = - Determinant.alfa a - |> Spindeterminant.bitstring - in - if Z.logand neg_active_mask alfa <> occ_mask then - false - else - let beta = - Determinant.beta a - |> Spindeterminant.bitstring - in - Z.logand neg_active_mask beta = occ_mask - in let symmetric = i_o1_alfa == alfa_o2_i in @@ -606,6 +580,35 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states } |> Util.stream_to_list + + +let is_internal det_space = + let m l = + List.fold_left (fun accu i -> + let j = i-1 in Z.(logor accu (shift_left one j)) + ) Z.zero l + in + let mo_class = DeterminantSpace.mo_class det_space in + let active_mask = m (MOClass.active_mos mo_class) in + let occ_mask = m (MOClass.core_mos mo_class) in + let inactive_mask = m (MOClass.inactive_mos mo_class) in + let occ_mask = Z.logor occ_mask inactive_mask in + let neg_active_mask = Z.lognot active_mask in + fun a -> + let alfa = + Determinant.alfa a + |> Spindeterminant.bitstring + in + if Z.logand neg_active_mask alfa <> occ_mask then + false + else + let beta = + Determinant.beta a + |> Spindeterminant.bitstring + in + Z.logand neg_active_mask beta = occ_mask + + let pt2_en ci = let mo_basis = Ds.mo_basis ci.det_space in @@ -652,7 +655,7 @@ let pt2_en ci = +. (fock_diag.(p) -. two_e h p h p s s) |> ignore; h_ij mo_basis alfa alfa - | _ -> assert false + | _ -> e0.{1} -. 1.0 in let e0 = e0.{1} in fun alfa -> @@ -667,12 +670,13 @@ let pt2_en ci = in second_order_sum ci list_holes list_particles list_holes list_particles - i_o1_alfa i_o1_alfa w_alfa psi0 + (is_internal ci.det_space) i_o1_alfa i_o1_alfa w_alfa psi0 |> List.fold_left (+.) 0. + let pt2_mp ci = let mo_basis = Ds.mo_basis ci.det_space in @@ -699,13 +703,14 @@ let pt2_mp ci = let psi0, _ = Parallel.broadcast ci.eigensystem in second_order_sum ci list_holes list_particles list_holes list_particles - i_o1_alfa i_o1_alfa w_alfa psi0 + (is_internal ci.det_space) i_o1_alfa i_o1_alfa w_alfa psi0 |> List.fold_left (+.) 0. let variance ci = let mo_basis = Ds.mo_basis ci.det_space in + let psi0, _ = Parallel.broadcast ci.eigensystem in let i_o1_alfa = h_ij mo_basis in @@ -718,9 +723,60 @@ let variance ci = [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ] in - let psi0, _ = Parallel.broadcast ci.eigensystem in second_order_sum ci list_holes list_particles list_holes list_particles - i_o1_alfa i_o1_alfa w_alfa psi0 + (is_internal ci.det_space) i_o1_alfa i_o1_alfa w_alfa psi0 |> List.fold_left (+.) 0. + + +let pt2_en_reference ci = + + let mo_basis = Ds.mo_basis ci.det_space in + let psi0, e0 = Parallel.broadcast ci.eigensystem in + + let aux_basis = mo_basis in + let ds = + DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis + in + let out_dets = + ds + |> DeterminantSpace.determinants_array + |> Array.to_list + |> List.filter (fun i -> not (is_internal ci.det_space i)) + |> Array.of_list + in + + let in_dets = + DeterminantSpace.determinants_array ci.det_space + in + + let m_H_aux = + let h_aa = + Array.map (fun ki -> h_ij aux_basis ki ki) out_dets + in + Array.map (fun ki -> + Array.mapi (fun j kj -> + (h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j)) + ) out_dets + ) in_dets + |> Mat.of_array + in + + let m_F_aux = + Array.map (fun ki -> + Array.map (fun kj -> + h_ij aux_basis ki kj + ) out_dets + ) in_dets + |> Mat.of_array + in + + let m_HF = + gemm m_H_aux m_F_aux ~transb:`T + in + (gemm ~transa:`T psi0 @@ gemm m_HF psi0).{1,1} + + + + diff --git a/CI/F12CI.ml b/CI/F12CI.ml index 061835d..52f1117 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -5,6 +5,7 @@ open Lacaml.D type t = { + gamma : float; mo_basis : MOBasis.t ; aux_basis : MOBasis.t ; det_space : DeterminantSpace.t ; @@ -20,14 +21,14 @@ let eigensystem t = Lazy.force t.eigensystem let f12_integrals mo_basis = - let two_e_ints = MOBasis.two_e_ints mo_basis in + let two_e_ints = MOBasis.f12_ints mo_basis in ( (fun i j _ -> 0.), (fun i j k l s s' -> if s' = Spin.other s then - ERI.get_phys two_e_ints i j k l + F12.get_phys two_e_ints i j k l else - (ERI.get_phys two_e_ints i j k l) -. - (ERI.get_phys two_e_ints i j l k) + (F12.get_phys two_e_ints i j k l) -. + (F12.get_phys two_e_ints i j l k) ) ) @@ -51,13 +52,36 @@ let f_ij mo_basis ki kj = |> List.hd -let dressing_vector f12_amplitudes ci = +let is_internal det_space = + let m l = + List.fold_left (fun accu i -> + let j = i-1 in Z.(logor accu (shift_left one j)) + ) Z.zero l + in + let mo_class = DeterminantSpace.mo_class det_space in + let aux_mask = m (MOClass.auxiliary_mos mo_class) in + fun a -> + let alfa = + Determinant.alfa a + |> Spindeterminant.bitstring + in + if Z.logand aux_mask alfa <> Z.zero then + false + else + let beta = + Determinant.beta a + |> Spindeterminant.bitstring + in + Z.logand aux_mask beta = Z.zero + + +let dressing_vector aux_basis 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 - let alfa_o2_i = f_ij mo_basis in +(* + let i_o1_alfa = h_ij aux_basis in + let alfa_o2_i = f_ij aux_basis in let w_alfa _ _ = 1. in @@ -65,19 +89,59 @@ debug "Computing dressing vector"; let list_holes = List.concat [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ] - and list_particles1 = MOClass.auxiliary_mos mo_class - and list_particles2 = List.concat + and list_particles2 = MOClass.auxiliary_mos mo_class + and list_particles1 = 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 ~unique:false + (is_internal ci.det_space) 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 |] +*) + + let out_dets = + DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis + |> DeterminantSpace.determinants_array + |> Array.to_list + |> List.filter (fun i -> not (is_internal ci.CI.det_space i)) + |> Array.of_list + in + + let in_dets = + DeterminantSpace.determinants_array ci.CI.det_space + in + + let m_H_aux = + Array.map (fun ki -> + Array.map (fun kj -> + h_ij aux_basis ki kj + ) out_dets + ) in_dets + |> Mat.of_array + in + + let m_F_aux = + Array.map (fun ki -> + Array.map (fun kj -> + f_ij aux_basis ki kj + ) out_dets + ) in_dets + |> Mat.of_array + in + +Printf.printf "%d %d\n" (Mat.dim1 m_H_aux) (Mat.dim2 m_H_aux); + let m_HF = + gemm m_H_aux m_F_aux ~transb:`T + in + gemm m_HF f12_amplitudes + |> Matrix.sparse_of_mat @@ -85,6 +149,8 @@ Util.debug_matrix "f12 amplitudes" f12_amplitudes; let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basis_filename () = + let gamma = 1. in + let mo_num = MOBasis.size mo_basis in (* Add auxiliary basis set *) @@ -118,6 +184,7 @@ let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basi Parallel.broadcast (lazy x) in + let f12_amplitudes = (* While in a sequential region, initiate the parallel 4-idx transformation to avoid nested parallel jobs @@ -127,9 +194,9 @@ debug "Four-idx transform of f12 intergals"; let f = fun ki kj -> if ki <> kj then - f_ij aux_basis ki kj + gamma *. (f_ij aux_basis ki kj) else - f_ij aux_basis ki kj +. 1. + 1. +. gamma *. (f_ij aux_basis ki kj) in debug "Computing F matrix"; let m_F = @@ -137,7 +204,10 @@ debug "Computing F matrix"; |> Lazy.force in fun ci_coef -> +(* +Util.debug_matrix "F" (Matrix.to_mat m_F); debug "Solving linear system"; +*) Matrix.ax_eq_b m_F (Matrix.dense_of_mat ci_coef) |> Matrix.to_mat in @@ -154,32 +224,43 @@ debug "Solving linear system"; let m_H = Lazy.force ci.CI.m_H in +(* +Util.debug_matrix "H" (Matrix.to_mat m_H); +*) let rec iteration ?(state=1) psi = +(* debug "Iteration"; +Util.debug_matrix "T" (f12_amplitudes psi); +*) let delta = - dressing_vector (f12_amplitudes psi) ci + dressing_vector aux_basis (f12_amplitudes psi) ci in 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) ; + let f = 1.0 /. psi.{1,1} in 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})); + |> List.iter (fun i -> + let delta = (Matrix.get delta i 1) *. f in + m_H_dressed.{i,1} <- m_H_dressed.{i,1} +. delta; + if i <> 1 then + begin + m_H_dressed.{1,i} <- m_H_dressed.{1,i} +. delta; + m_H_dressed.{1,1} <- m_H_dressed.{1,1} -. delta *. psi.{i,1} *. f; + end + ); 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)) + 1.0 -. abs_float ( dot + (Mat.to_col_vecs psi).(0) + (Mat.to_col_vecs eigenvectors).(0) ) in - Printf.printf "Convergence : %f %f\n" conv eigenvalues.{1}; + Printf.printf "Convergence : %f %f\n" conv (eigenvalues.{1} +. e_shift); if conv > threshold then iteration eigenvectors else @@ -193,6 +274,7 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta; ------- *) (* +(*TODO : enlever le double comptage de la symmetrisation*) 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.) ) @@ -226,6 +308,6 @@ Format.printf "Dressing vector : %a\n@." Matrix.pp_matrix delta; ) in - { mo_basis ; aux_basis ; det_space ; ci ; eigensystem } + { mo_basis ; aux_basis ; det_space ; ci ; eigensystem ; gamma } diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index 376d511..490e3ae 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -303,6 +303,11 @@ let four_index_transform coef source = incr jk; get_chem_all_i source ~j ~k ~l |> Array.iteri (fun i x -> o.{i+1,!jk} <- x) + (* + lacpy ~bc:!jk ~b:o + (Mat.of_col_vecs [| Vec.of_array (get_chem_all_i source ~j ~k ~l) |] ) + |> ignore + *) ) range_ao ) range_ao; (* o_i_jk *) diff --git a/run_cas.ml b/run_cas.ml index 75266ea..a58a330 100644 --- a/run_cas.ml +++ b/run_cas.ml @@ -72,14 +72,22 @@ let () = let ci = CI.make space in Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); +(* let pt2 = CI.pt2_mp ci in Format.fprintf ppf "CAS-MP2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2); +*) let pt2 = CI.pt2_en ci in - Format.fprintf ppf "CAS-EN2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2); + Format.fprintf ppf "CAS-EN2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2); + let pt2 = CI.pt2_en_reference ci in + Format.fprintf ppf "CAS-EN2 energy (reference) : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2); + +(* let variance = CI.variance ci in Format.fprintf ppf "CAS variance : %20.16f@." variance + *) + (* let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in Util.list_range 1 (DeterminantSpace.size space)