diff --git a/CI/F12CI.ml b/CI/F12CI.ml index b4f6704..54f1356 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -20,11 +20,16 @@ let f12_integrals mo_basis = let two_e_ints = MOBasis.f12_ints mo_basis in ( (fun _ _ _ -> 0.), (fun i j k l s s' -> - if s' = Spin.other s then - F12.get_phys two_e_ints i j k l + if (i=k && j<>l) || (j=l && i<>k) then + 0. else - (F12.get_phys two_e_ints i j k l) -. - (F12.get_phys two_e_ints i j l k) + begin + if s' = Spin.other s then + 0.5 *. (F12.get_phys two_e_ints i j k l) + else + 0.25 *. ((F12.get_phys two_e_ints i j k l) -. + (F12.get_phys two_e_ints i j l k)) + end ) ) @@ -70,16 +75,18 @@ let is_internal det_space = Determinant.alfa a |> Spindeterminant.bitstring in - (* let beta = Determinant.beta a |> Spindeterminant.bitstring in let a = Bitstring.logand aux_mask alfa and b = Bitstring.logand aux_mask beta - in Bitstring.popcount a + Bitstring.popcount b < 2 - *) + in + match Bitstring.popcount a + Bitstring.popcount b with + | 1 | 2 -> false + | _ -> true + (* if not (Bitstring.logand aux_mask alfa |> Bitstring.is_zero ) then false else @@ -88,6 +95,7 @@ let is_internal det_space = |> Spindeterminant.bitstring in Bitstring.logand aux_mask beta |> Bitstring.is_zero + *) let dressing_vector ~frozen_core aux_basis f12_amplitudes ci = @@ -118,49 +126,88 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci = Printf.printf "Building matrix\n%!"; + (* Determinants of the FCI space as a list *) + let in_dets = + DeterminantSpace.determinants_array ci.CI.det_space + |> Array.to_list + in - let m_H_aux, m_F_aux = - let out_dets_stream = + (* Stream that generates only singly and doubly excited determinants + wrt FCI space *) + let out_dets_stream = + + (* Stream that generates all determinants of FCI space *) + let s = DeterminantSpace.fci_of_mo_basis ~frozen_core aux_basis |> DeterminantSpace.determinant_stream - in - let in_dets = - DeterminantSpace.determinants_array ci.CI.det_space - |> Array.to_list in + (* Select only singly and doubly excited determinants + wrt FCI space *) + Stream.from (fun _ -> + try + let rec result () = + let ki = Stream.next s in + if is_internal ci.CI.det_space ki then + result () + else + Some ki + in + result () + with Stream.Failure -> None + ) + in - let rec col_vecs_list accu_H accu_F = + let make_h_and_f n = + + let rec col_vecs_list accu_H accu_F = function + | 0 -> + List.rev accu_H, + List.rev accu_F + | n -> try let ki = Stream.next out_dets_stream in - if is_internal ci.CI.det_space ki then - raise Exit - else - let h, f = - List.map (fun kj -> - match hf_ij aux_basis ki kj with - | [ a ; b ] -> a, b - | _ -> assert false ) in_dets - |> List.split - in - col_vecs_list (h::accu_H) (f::accu_F) + let h, f = + List.map (fun kj -> + match hf_ij aux_basis ki kj with + | [ a ; b ] -> a, b + | _ -> assert false ) in_dets + |> List.split + in + let h = + Vec.of_list h + and f = + Vec.of_list f + in + col_vecs_list (h::accu_H) (f::accu_F) (n-1) with - | Exit -> col_vecs_list accu_H accu_F - | Stream.Failure -> - List.rev_map Vec.of_list accu_H, - List.rev_map Vec.of_list accu_F + | Stream.Failure -> col_vecs_list accu_H accu_F 0 in let h, f = - col_vecs_list [] [] + col_vecs_list [] [] n in - Mat.of_col_vecs_list h, + Mat.of_col_vecs_list h, Mat.of_col_vecs_list f in Printf.printf "Matrix product\n%!"; let m_HF = - gemm m_H_aux m_F_aux ~transb:`T + let batch_size = 10_000_000 / (Mat.dim1 f12_amplitudes) in + let result = + let m_H_aux, m_F_aux = make_h_and_f batch_size in + gemm m_H_aux m_F_aux ~transb:`T + in + while (Stream.peek out_dets_stream <> None) + do + Printf.printf "gemm\n%!"; + let m_H_aux, m_F_aux = make_h_and_f batch_size in + let hf = + gemm m_H_aux m_F_aux ~transb:`T + in + ignore @@ Mat.add result hf ~c:result + done; + result in @@ -200,7 +247,7 @@ let dressing_vector ~frozen_core aux_basis f12_amplitudes ci = Printf.printf "Done\n%!"; gemm m_HF f12_amplitudes - |> Matrix.sparse_of_mat + |> Matrix.dense_of_mat @@ -211,6 +258,7 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen let f12 = Util.of_some @@ Simulation.f12 simulation in let mo_num = MOBasis.size mo_basis in +Printf.printf "Add aux basis\n%!"; (* Add auxiliary basis set *) let s = let charge = Charge.to_int @@ Simulation.charge simulation @@ -230,7 +278,16 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen let aux_basis = MOBasis.of_mo_basis s mo_basis in + let () = +Printf.printf "F12 ints\n%!"; + ignore @@ MOBasis.f12_ints aux_basis + in + let () = +Printf.printf "2e ints\n%!"; + ignore @@ MOBasis.two_e_ints aux_basis + in +Printf.printf "det space\n%!"; let det_space = DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num in @@ -243,12 +300,12 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen in + (* let f12_amplitudes = fun c -> let result = lacpy c in Mat.scal (0.5) result ; result - (* (* While in a sequential region, initiate the parallel 4-idx transformation to avoid nested parallel jobs *) @@ -276,8 +333,8 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen if abs_float norm.{1,1} > 1. then failwith "Norm of |F> > 1"; result - *) in + *) let e_shift = let det = @@ -294,7 +351,7 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen let rec iteration ?(state=1) psi = let delta = - dressing_vector ~frozen_core aux_basis (f12_amplitudes psi) ci + dressing_vector ~frozen_core aux_basis psi ci in let f = 1.0 /. psi.{1,1} in diff --git a/run_fci_f12.ml b/run_fci_f12.ml index c65a598..5a374a0 100644 --- a/run_fci_f12.ml +++ b/run_fci_f12.ml @@ -76,7 +76,7 @@ let () = in let fcif12 = - F12CI.make ~simulation ~frozen_core:true ~mo_basis ~aux_basis_filename () + F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename () in let ci = F12CI.ci fcif12 in