(** %{ $ \langle ij | H F | kl \rangle $ %} integrals. *)

open Lacaml.D

module Fis = FourIdxStorage


type t = {
  simulation : Simulation.t ;
  aux_basis  : MOBasis.t ;
  f_0 : Determinant.t -> float ;
  f_1 : Determinant.t -> Determinant.t -> float ;
  f_2 : Determinant.t -> Determinant.t -> float ;
  f_3 : Determinant.t -> Determinant.t -> float ;
}


let sum l f = Array.fold_left (fun accu i -> accu +. f i) 0. l

let array_3_init d1 d2 d3 fx =
  let f k = 
    let x = 
      Bigarray.(Array2.create Float64 fortran_layout d1 d2)
    in
    for j=1 to d2 do
      for i=1 to d1 do
        x.{i,j} <-  fx i j k
      done
    done;
    (k,x)
  in
  let result =
      SharedMemory.create Bigarray.Float64 [| d1 ; d2 ; d3 |]
      |> Bigarray.array3_of_genarray
  in
  Util.list_range 1 d3
  |> Stream.of_list 
  |> Farm.run ~f ~ordered:false
  |> Stream.iter (fun (k,x) ->
        for j=1 to d2 do
            for i=1 to d1 do
              result.{i,j,k} <-  x.{i,j}
            done
        done)
  ;
  if Parallel.master then Printf.printf "Broadcast d3\n%!" ;
  try Parallel.InterNode.broadcast (lazy result) 
  with Invalid_argument _ ->
   begin
      Printf.eprintf "Array too large.\n%!";
      failwith "Array too large"
   end

let array_4_init d1 d2 d3 d4 fx =
  let f (k,l) = 
    let x = 
      Bigarray.(Array2.create Float64 fortran_layout d1 d2)
    in
    for j=1 to d2 do
      for i=1 to d1 do
        x.{i,j} <-  fx i j k l
      done
    done;
    (k,l,x)
  in
  let result =
      SharedMemory.create Bigarray.Float64 [| d1;d2;d3;d4 |] 
  in
  Util.list_range 1 d4
  |> List.rev_map (fun l -> 
    Util.list_range 1 d3
    |> List.rev_map (fun k -> (k,l)) )
  |> List.concat
  |> List.rev
  |> Stream.of_list 
  |> Farm.run ~f ~ordered:false
  |> Stream.iter (fun (k,l,x) ->
        for j=1 to d2 do
          for i=1 to d1 do
            result.{i,j,k,l} <-  x.{i,j}
          done
        done)
  ;
  if Parallel.master then Printf.printf "Broadcast d4\n%!" ;
  try Parallel.InterNode.broadcast (lazy result)
  with Invalid_argument _ ->
    begin
      Printf.eprintf "Array too large... splitting.\n%!";
      let x = 
        Bigarray.(Array3.create Float64 fortran_layout d1 d2 d3)
      in
      for l=1 to d4 do
        if Parallel.master then
          begin
            for k=1 to d3 do
              for j=1 to d2 do
                for i=1 to d1 do
                  x.{i,j,k} <- result.{i,j,k,l} 
                done
              done
            done;
            ignore @@ Parallel.InterNode.broadcast (lazy x)
          end
        else 
          begin
            ignore @@ Parallel.InterNode.broadcast (lazy x);
            for k=1 to d3 do
              for j=1 to d2 do
                for i=1 to d1 do
                  result.{i,j,k,l} <- x.{i,j,k} 
                done
              done
            done
          end
      done;
      result
    end

let array_5_init d1 d2 d3 d4 d5 fx =
  let f (l,m) =
    let x = 
      Bigarray.(Array3.create Float64 fortran_layout d1 d2 d3)
    in
    for k=1 to d3 do
      for j=1 to d2 do
        for i=1 to d1 do
          x.{i,j,k} <- fx i j k l m
        done
      done
    done;
    (l,m,x)
  in
  let result =
      SharedMemory.create Bigarray.Float64 [| d1;d2;d3;d4;d5 |]
  in
  Util.list_range 1 d5
  |> List.rev_map (fun m -> 
    Util.list_range 1 d4
    |> List.rev_map (fun l -> (l,m)) )
  |> List.concat
  |> List.rev
  |> Stream.of_list 
  |> Farm.run ~f ~ordered:false
  |> Stream.iter (fun (l,m,x) ->
      for k=1 to d3 do
        for j=1 to d2 do
          for i=1 to d1 do
            result.{i,j,k,l,m} <-  x.{i,j,k}
          done
        done
      done)
  ;
  if Parallel.master then Printf.printf "Broadcast d5\n%!" ;
  try
    Parallel.InterNode.broadcast (lazy result)
  with Invalid_argument _ ->
    begin
      Printf.eprintf "Array too large... splitting.\n%!";
      let x = 
        Bigarray.(Genarray.create Float64 fortran_layout [| d1; d2; d3; d4 |])
      in
      for m=1 to d5 do
        if Parallel.master then
          begin
            for l=1 to d4 do
              for k=1 to d3 do
                for j=1 to d2 do
                  for i=1 to d1 do
                    x.{i,j,k,l} <- result.{i,j,k,l,m} 
                  done
                done
              done
            done;
            ignore @@ Parallel.InterNode.broadcast (lazy x)
          end
        else 
          begin
            ignore @@ Parallel.InterNode.broadcast (lazy x);
            for l=1 to d4 do
              for k=1 to d3 do
                for j=1 to d2 do
                  for i=1 to d1 do
                    result.{i,j,k,l,m} <- x.{i,j,k,l} 
                  done
                done
              done
            done
          end
      done;
      result
    end


let make ~frozen_core ~simulation ~mo_basis ~aux_basis_filename () =
  
  let f12 = Util.of_some @@ Simulation.f12 simulation in                                      
  let mo_num = MOBasis.size mo_basis in

  (* Add auxiliary basis set *)                                                               
  let simulation =
    let charge        = Charge.to_int @@ Simulation.charge simulation 
    and multiplicity  = Electrons.multiplicity @@ Simulation.electrons simulation
    and nuclei        = Simulation.nuclei simulation
    in
    let general_basis = 
      Basis.general_basis @@ Simulation.basis simulation
    in
    GeneralBasis.combine [
      general_basis ; GeneralBasis.read aux_basis_filename
    ]
    |> Basis.of_nuclei_and_general_basis ~f12 nuclei
    |> Simulation.make ~charge ~multiplicity ~nuclei 
  in

  let aux_basis = 
    MOBasis.of_mo_basis simulation mo_basis
  in

  let aux_num = MOBasis.size aux_basis in

  (* While in a sequential region, initiate the parallel
     4-idx transformation to avoid nested parallel jobs
  *)
  ignore @@ MOBasis.two_e_ints aux_basis;
  ignore @@ MOBasis.f12_ints aux_basis;

  (* Compute the <ij|QHF|kl> integrals *)
  if Parallel.master then Printf.eprintf "Computing HF12 integrals\n%!";

  let mos_cabs =
    Util.list_range (mo_num+1) aux_num
    |> Array.of_list 
  in

  let n_core =
    if frozen_core then
      (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
    else 0
  in

  let mos_in   =
    Util.list_range (n_core+1) mo_num
    |> Array.of_list
  in

  let mos_a k =
    Determinant.alfa k
    |> Spindeterminant.to_list
    |> List.filter (fun i -> i > n_core)
    |> Array.of_list 
  in

  let mos_b k =
    Determinant.beta k
    |> Spindeterminant.to_list
    |> List.filter (fun i -> i > n_core)
    |> Array.of_list 
  in

  let h_one =
    let h =
      MOBasis.one_e_ints aux_basis
    in fun i j _ -> h.{i,j}
  in

  let h_two =
    let two_e_ints = MOBasis.two_e_ints aux_basis in
    let h2 i j k l (s:Spin.t) (s':Spin.t) =
      if s' <> s then
          ERI.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)
    in
    h2
  in

  let f_two =
    let f12_ints = MOBasis.f12_ints aux_basis in
    let f2 i j k l (s:Spin.t) (s':Spin.t) =
      if s' <> s then
          0.375 *. F12.get_phys f12_ints i j k l +.
          0.125 *. F12.get_phys f12_ints i j l k
      else
          0.25 *. (
          (F12.get_phys f12_ints i j k l) -.
          (F12.get_phys f12_ints i j l k) )
    in
    f2
  in

  let f_one = fun _ _ _ -> 0. in


(*
  let s, s' = Spin.(Alfa, Beta) in
  let m_3_Hschwarz_aa = 
    array_3_init mo_num mo_num mo_num (fun i j k ->
      sum mos_cabs (fun a -> let x = h_two a i j k s s in x *. x)
    )
  in
  let m_3_Hschwarz_ab = 
    array_3_init mo_num mo_num mo_num (fun i j k ->
      sum mos_cabs (fun a -> let x = h_two a i j k s s' in x *. x)
    )
  in
  let m_3_Fschwarz_aa = 
    array_3_init mo_num mo_num mo_num (fun i j k ->
      sum mos_cabs (fun a -> let x = f_two a i j k s s in x *. x)
    )
  in
  let m_3_Fschwarz_ab = 
    array_3_init mo_num mo_num mo_num (fun i j k ->
      sum mos_cabs (fun a -> let x = f_two a i j k s s' in x *. x)
    )
  in

  let eps = 1.e-15 in

  let sum_schwarz i j k l m n s1 s2 s3 s4 =
    if (s1  = s2 && m_3_Hschwarz_aa.{i,j,k} < eps ) ||
       (s1 <> s2 && m_3_Hschwarz_ab.{i,j,k} < eps ) ||
       (s3  = s4 && m_3_Fschwarz_aa.{l,m,n} < eps ) ||
       (s3 <> s4 && m_3_Fschwarz_ab.{l,m,n} < eps ) 
    then
      0.
    else
      Array.fold_left (fun accu a ->
          accu +. h_two a i j k s1 s2 *. f_two a l m n s3 s4
          ) 0. mos_cabs
  in
*)

  (* Pre-compute dressed integrals *)
  if Parallel.master then Printf.printf "Computing m_0111_1H_1F\n%!" ;
  let m_0111_1H_1F =
    Vec.init mo_num (fun i ->
          sum mos_cabs (fun a ->
              h_one a i Spin.Alfa *. f_one a i Spin.Alfa ))
  in

  if Parallel.master then Printf.printf "Computing m_0111_1H_2Fa\n%!" ;
  let m_0111_1H_2Fa, m_0111_2Ha_2Fa =

      let m_0122_Haa =
        array_3_init mo_num mo_num mo_num (fun i j k ->
            sum mos_cabs (fun a ->
              h_two a k i j Spin.Alfa Spin.Alfa *. f_two a k i j Spin.Alfa Spin.Alfa
              ) )
      in

      let m_0111_1H_2Fa =
        Mat.init_cols mo_num mo_num (fun i j ->
            sum mos_cabs (fun a ->
              h_one a i Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Alfa +.
              h_two a j i j Spin.Alfa Spin.Alfa *. f_one a i Spin.Alfa
            ) +.
            if i < j then 0. else
                begin
                    sum mos_cabs (fun a ->
                        sum mos_cabs (fun b -> if b >= a then 0. else
                          h_two i j a b Spin.Alfa Spin.Alfa *. f_two a b i j Spin.Alfa Spin.Alfa
                        )
                    ) +.
                    sum mos_in (fun k -> m_0122_Haa.{i,j,k})
                end
        )
      in

      let m_0111_2Ha_2Fa =
        array_3_init mo_num mo_num mo_num (fun i j k ->
          sum mos_cabs (fun a ->
            h_two a j i j Spin.Alfa Spin.Alfa *.
            f_two a k i k Spin.Alfa Spin.Alfa
          ) -. if i < j then 0. else m_0122_Haa.{i,j,k}
        )
      in m_0111_1H_2Fa, m_0111_2Ha_2Fa
  in

  if Parallel.master then Printf.printf "Computing m_0111_1H_2Fb\n%!" ;
  let m_0111_1H_2Fb, m_0111_2Hb_2Fb =
      let m_0122_Hab =
        array_3_init mo_num mo_num mo_num (fun i j k ->
            sum mos_cabs (fun a ->
              h_two a k i j Spin.Alfa Spin.Beta *. f_two a k i j Spin.Alfa Spin.Beta
            ) )
    in

    let m_0111_1H_2Fb =
      Mat.init_cols mo_num mo_num (fun i j ->
          sum mos_cabs (fun a ->
            h_one a i Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Beta +.
            h_two a j i j Spin.Alfa Spin.Beta *. f_one a i Spin.Alfa +.
            h_one a j Spin.Alfa *. f_two a i j i Spin.Alfa Spin.Beta +.
            h_two a i j i Spin.Alfa Spin.Beta *. f_one a j Spin.Alfa
          ) +.
        sum mos_in (fun k -> m_0122_Hab.{i,j,k} +. m_0122_Hab.{j,i,k} ) +.
        sum mos_cabs (fun b ->
          sum mos_cabs (fun a ->
            h_two a b i j Spin.Alfa Spin.Beta *. f_two a b i j Spin.Alfa Spin.Beta
            )
          )
        )
    in

    let m_0111_2Hb_2Fb =
      array_3_init mo_num mo_num mo_num (fun i j k ->
        sum mos_cabs (fun a ->
          h_two a i k i Spin.Alfa Spin.Beta *.
          f_two a j k j Spin.Alfa Spin.Beta +.
          h_two a k i k Spin.Alfa Spin.Beta *.
          f_two a j i j Spin.Alfa Spin.Alfa
        ) -. m_0122_Hab.{k,i,j}
      )
    in
    m_0111_1H_2Fb, m_0111_2Hb_2Fb
  in

  if Parallel.master then Printf.printf "Computing m_0111_2Ha_2Fb\n%!" ;
  let m_0111_2Ha_2Fb =
    array_3_init mo_num mo_num mo_num (fun i j k ->
      sum mos_cabs (fun a ->
        h_two a j i j Spin.Alfa Spin.Alfa *.
        f_two a k i k Spin.Alfa Spin.Beta
      )
    )
  in

  let f_0 ki =

    let mos_i, mos_i' = mos_a ki, mos_b ki in
    let same = (mos_i = mos_i') in

    (* Alpha *)
    let a =
        sum mos_i  (fun i -> m_0111_1H_1F.{i})
    in

    let b =
        if same then a else
        sum mos_i' (fun i -> m_0111_1H_1F.{i})
    in

    let aa =
        sum mos_i (fun j ->
          sum mos_i (fun i -> m_0111_1H_2Fa.{i,j} ))
    in

    let bb =
      if same then aa else
      sum mos_i' (fun j ->
        sum mos_i' (fun i -> m_0111_1H_2Fa.{i,j} ))
    in

    let ab =
      sum mos_i' (fun j ->
        sum mos_i (fun i -> m_0111_1H_2Fb.{i,j}  ))
    in

    let aaa =
        sum mos_i (fun k ->
          sum mos_i (fun j ->
            sum mos_i (fun i -> m_0111_2Ha_2Fa.{i,j,k} )))
    in

    let bbb =
      if same then aaa else
        sum mos_i' (fun k ->
          sum mos_i' (fun j ->
            sum mos_i' (fun i -> m_0111_2Ha_2Fa.{i,j,k} )))
    in

    let baa =
        sum mos_i' (fun k ->
          sum mos_i (fun j ->
            sum mos_i (fun i ->
              m_0111_2Ha_2Fb.{i,j,k} +. m_0111_2Hb_2Fb.{i,j,k}
        )))
    in

    let bba =
        sum mos_i  (fun k ->
          sum mos_i' (fun j ->
            sum mos_i' (fun i ->
              m_0111_2Ha_2Fb.{i,j,k} +. m_0111_2Hb_2Fb.{j,i,k}
        )))
    in

    a +. b +. aa +. bb +. ab +. aaa +. baa +. bba +. bbb
  in

  if Parallel.master then Printf.printf "Computing m_1111_1H_1F\n%!" ;
  let m_1111_1H_1F =
    Mat.init_cols mo_num mo_num (fun i k ->
          sum mos_cabs (fun a -> h_one a i Spin.Alfa *. f_one a k Spin.Alfa ))
  in

  if Parallel.master then Printf.printf "Computing m_1111_2Ha_2Fa\n%!" ;
  let m_1111_2Ha_2Fa =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      if l=i then
          sum mos_cabs (fun a ->
            h_two a l j l Spin.Alfa Spin.Alfa *.
            f_two a i j k Spin.Alfa Spin.Alfa )
      else
          sum mos_cabs (fun a ->
            h_two a j i j Spin.Alfa Spin.Alfa *.
            f_two a l k l Spin.Alfa Spin.Alfa +.
            h_two a l j l Spin.Alfa Spin.Alfa *.
            f_two a i j k Spin.Alfa Spin.Alfa )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1111_2Hb_2Fa\n%!" ;
  let m_1111_2Hb_2Fa =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      if l=i then
          sum mos_cabs (fun a ->
            h_two a l j l Spin.Alfa Spin.Beta *.
            f_two a i j k Spin.Alfa Spin.Beta )
      else
          sum mos_cabs (fun a ->
            h_two a j i j Spin.Alfa Spin.Beta *.
            f_two a l k l Spin.Alfa Spin.Alfa +.
            h_two a l j l Spin.Alfa Spin.Beta *.
            f_two a i j k Spin.Alfa Spin.Beta )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1111_2Ha_2Fb\n%!" ;
  let m_1111_2Ha_2Fb =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a j i j Spin.Alfa Spin.Alfa *.
        f_two a l k l Spin.Alfa Spin.Beta +.
        h_two a l j l Spin.Alfa Spin.Beta *.
        f_two a i j k Spin.Alfa Spin.Alfa 
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1111_2Hb_2Fb\n%!" ;
  let m_1111_2Hb_2Fb =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a j i j Spin.Alfa Spin.Beta *.
        f_two a l k l Spin.Alfa Spin.Beta +.
        h_two a l j l Spin.Alfa Spin.Alfa *.
        f_two a i j k Spin.Alfa Spin.Beta 
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1121_2Ha_2Fa\n%!" ;
  let m_1121_2Ha_2Fa =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a k j i Spin.Alfa Spin.Alfa *.
        f_two a l j l Spin.Alfa Spin.Alfa 
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1121_2Hb_2Fa\n%!" ;
  let m_1121_2Hb_2Fa =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a k j i Spin.Alfa Spin.Beta *.
        f_two a l j l Spin.Alfa Spin.Alfa 
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1121_2Ha_2Fb\n%!" ;
  let m_1121_2Ha_2Fb =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a k j i Spin.Alfa Spin.Alfa *.
        f_two a l j l Spin.Alfa Spin.Beta 
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1121_2Hb_2Fb\n%!" ;
  let m_1121_2Hb_2Fb =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a k j i Spin.Alfa Spin.Beta *.
        f_two a l j l Spin.Alfa Spin.Beta 
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1122_va\n%!" ;
  let m_1122_va =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a l j i Spin.Alfa Spin.Alfa *.
        f_two a l j k Spin.Alfa Spin.Alfa 
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1111_1H_2Fa\n%!" ;
  let m_1111_1H_2Fa =
    array_3_init mo_num mo_num mo_num (fun j i k ->
        sum mos_in (fun l -> m_1122_va.{l,j,i,k}) +.
        sum mos_cabs (fun a ->
          h_one a i Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Alfa +.
          h_two a j i j Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.
          h_one a j Spin.Alfa *. f_two a i j k Spin.Alfa Spin.Alfa +.
          h_two a k j i Spin.Alfa Spin.Alfa *. f_one a j Spin.Alfa +.
          sum mos_cabs (fun b -> if b > a then 0. else
            h_two b a j i Spin.Alfa Spin.Alfa *.
            f_two b a j k Spin.Alfa Spin.Alfa )
        )
      )
  in

  if Parallel.master then Printf.printf "Computing m_1122_v2\n%!" ;
  let m_1122_v2 =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a l j i Spin.Alfa Spin.Beta *.
        f_two a l j k Spin.Alfa Spin.Beta 
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1122_v3\n%!" ;
  let m_1122_v3 =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a l i j Spin.Alfa Spin.Beta *.
        f_two a l k j Spin.Alfa Spin.Beta
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_1111_1H_2Fb\n%!" ;
  let m_1111_1H_2Fb =
    array_3_init mo_num mo_num mo_num (fun j i k ->
      sum mos_in (fun l -> m_1122_v2.{l,j,i,k} +. m_1122_v3.{l,j,i,k}) +.
      sum mos_cabs (fun a ->
        h_one a i Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Beta +.
        h_two a j i j Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.
        h_one a j Spin.Beta *. f_two a i j k Spin.Alfa Spin.Beta +.
        h_two a k j i Spin.Alfa Spin.Beta *. f_one a j Spin.Beta +.
        sum mos_cabs (fun b ->
          h_two b a j i Spin.Alfa Spin.Beta *.
          f_two b a j k Spin.Alfa Spin.Beta 
        )
      )
      )
  in

  if Parallel.master then Printf.printf "Computing m_1122_oa\n%!" ;
  let m_1122_oa =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
    if l > j then
      sum mos_cabs (fun a ->
        h_two a k j l Spin.Alfa Spin.Alfa *.
        f_two a i j l Spin.Alfa Spin.Alfa 
      )
    else 0.
    )
  in

  if Parallel.master then Printf.printf "Computing m_1122_o\n%!" ;
  let m_1122_o =
    array_4_init mo_num mo_num mo_num mo_num (fun l j i k ->
      sum mos_cabs (fun a ->
        h_two a k j l Spin.Alfa Spin.Beta *.
        f_two a i j l Spin.Alfa Spin.Beta 
      )
    )
  in

  let f_1 ki kj =
      let i, k, s, phase =
            match Excitation.of_det ki kj with
            | Excitation.(Single (phase, { hole ; particle ; spin })) ->
              hole, particle, spin, phase
            | _ -> assert false
      in

      let mos_novirt, mos_novirt' =
        let alfa =
          let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
          let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
          Bitstring.to_list (Bitstring.logor i j)
          |> Array.of_list 
        in
        let beta  =
          let i = Spindeterminant.bitstring @@ Determinant.beta ki in
          let j = Spindeterminant.bitstring @@ Determinant.beta kj in
          Bitstring.to_list (Bitstring.logor i j)
          |> Array.of_list 
        in
        match s with
        | Spin.Alfa -> alfa, beta
        | Spin.Beta -> beta, alfa
      in

      let mos_ij, mos_ij' =
        let alfa =
          let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
          let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
          Bitstring.to_list (Bitstring.logand i j)
          |> Array.of_list 
        in
        let beta  =
          let i = Spindeterminant.bitstring @@ Determinant.beta ki in
          let j = Spindeterminant.bitstring @@ Determinant.beta kj in
          Bitstring.to_list (Bitstring.logand i j)
          |> Array.of_list 
        in
        match s with
        | Spin.Alfa -> alfa, beta
        | Spin.Beta -> beta, alfa
      in

      let mos_i, mos_i' =
        match s with
        | Spin.Alfa -> mos_a ki, mos_b ki
        | Spin.Beta -> mos_b ki, mos_a ki
      in

      let mos_j, mos_j' =
        match s with
        | Spin.Alfa -> mos_a kj, mos_b kj
        | Spin.Beta -> mos_b kj, mos_a kj
      in

      let result =
      m_1111_1H_1F.{i,k} +.
      sum mos_ij (fun j ->
            m_1111_1H_2Fa.{j,i,k}
        +. sum mos_i  (fun l -> m_1111_2Ha_2Fa.{l,j,i,k})
        +. sum mos_i' (fun l -> m_1111_2Ha_2Fb.{l,j,i,k})
        +. sum mos_j  (fun l -> m_1121_2Ha_2Fa.{l,j,i,k})
        +. sum mos_j' (fun l -> m_1121_2Ha_2Fb.{l,j,i,k})
        -. sum mos_novirt (fun l -> m_1122_va.{l,j,i,k})
        -. sum mos_ij (fun l -> m_1122_oa.{l,j,i,k} )
        ) +.
      sum mos_ij' (fun j ->
            m_1111_1H_2Fb.{j,i,k}
        +. sum mos_i  (fun l -> m_1111_2Hb_2Fa.{l,j,i,k})
        +. sum mos_i' (fun l -> m_1111_2Hb_2Fb.{l,j,i,k})
        +. sum mos_j  (fun l -> m_1121_2Hb_2Fb.{l,j,i,k})
        +. sum mos_j' (fun l -> m_1121_2Hb_2Fa.{l,j,i,k})
        -. sum mos_novirt  (fun l ->  m_1122_v2.{l,j,i,k})
        -. sum mos_novirt' (fun l ->  m_1122_v3.{l,j,i,k})
        -. sum mos_ij    (fun l ->  m_1122_o.{l,j,i,k})
        )
      in
      match phase with
      | Phase.Pos -> result
      | Phase.Neg -> -. result
  in

  if Parallel.master then Printf.printf "Computing m_2112_1H_2Fa\n%!" ;
  let m_2112_1H_2Fa =
    array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
      sum mos_cabs (fun a ->
          h_one a i Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +.
          h_one a j Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Alfa +. 
          h_two a l i j Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.
          h_two a k j i Spin.Alfa Spin.Alfa *. f_one a l Spin.Alfa ) +.
      sum mos_cabs (fun a ->
          sum mos_in  (fun m -> -. h_two m a j i Spin.Alfa Spin.Alfa *.
                                  f_two m a k l Spin.Alfa Spin.Alfa) ) +.
      sum mos_cabs (fun a ->
          sum mos_cabs (fun b -> if b >= a then 0. else
              h_two b a j i Spin.Alfa Spin.Alfa *. f_two b a l k Spin.Alfa Spin.Alfa 
          )
        )
      )
  in

  if Parallel.master then Printf.printf "Computing m_2112_1H_2Fb\n%!" ;
  let m_2112_1H_2Fb =
    array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
      sum mos_cabs (fun a ->
          h_one a i Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.
          h_one a j Spin.Alfa *. f_two a i l k Spin.Alfa Spin.Beta +.
          h_two a l i j Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.
          h_two a k j i Spin.Alfa Spin.Beta *. f_one a l Spin.Alfa ) +.
      sum mos_cabs (fun a ->
          sum mos_in (fun m ->
            h_two m a j i Spin.Alfa Spin.Beta *.  f_two m a l k Spin.Alfa Spin.Beta +.
            h_two m a i j Spin.Alfa Spin.Beta *.  f_two m a k l Spin.Alfa Spin.Beta ) ) +.
      sum mos_cabs (fun a ->
          sum mos_cabs (fun b -> 
              h_two b a j i Spin.Alfa Spin.Beta *. f_two b a l k Spin.Alfa Spin.Beta 
          )
        )
      )
  in

  if Parallel.master then Printf.printf "Computing m_2112_2Ha_2Fa\n%!" ;
  let m_2112_2Ha_2Fa =
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a ->
          h_two a n i n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +. 
          h_two a n j n Spin.Alfa Spin.Alfa *. f_two a i l k Spin.Alfa Spin.Alfa 
        )
      )
  in

  if Parallel.master then Printf.printf "Computing m_2112_2Hb_2Fa\n%!" ;
  let m_2112_2Hb_2Fa =
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
    sum mos_cabs (fun a ->
        h_two a n i n Spin.Alfa Spin.Beta *. f_two a j k l Spin.Alfa Spin.Alfa +.
        h_two a n j n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Alfa
        )
      )
  in

  if Parallel.master then Printf.printf "Computing m_2112_2Ha_2Fb\n%!" ;
  let m_2112_2Ha_2Fb =
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
    sum mos_cabs (fun a ->
        h_two a n i n Spin.Alfa Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.
        h_two a n j n Spin.Alfa Spin.Beta *. f_two a i l k Spin.Alfa Spin.Beta
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_2121_2Ha_2Fa\n%!" ;
  let m_2121_2Ha_2Fa =
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a ->
          h_two a l i j Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Alfa +.
          h_two a k j i Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Alfa 
        )
      )
  in

  if Parallel.master then Printf.printf "Computing m_2121_2Hb_2Fa\n%!" ;
  let m_2121_2Hb_2Fa =
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a ->
          h_two a l i j Spin.Alfa Spin.Beta *. f_two a n k n Spin.Alfa Spin.Alfa +.
          h_two a k j i Spin.Alfa Spin.Beta *. f_two a n l n Spin.Alfa Spin.Beta 
        )
      )
  in

  if Parallel.master then Printf.printf "Computing m_2121_2Ha_2Fb\n%!" ;
  let m_2121_2Ha_2Fb =
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a ->
          h_two a l i j Spin.Alfa Spin.Alfa *. f_two a n k n Spin.Alfa Spin.Beta +.
          h_two a k j i Spin.Alfa Spin.Alfa *. f_two a n l n Spin.Alfa Spin.Beta 
        )
      )
  in

  if Parallel.master then Printf.printf "Computing m_2122_2Ha_2Fa_ij\n%!" ;
  let m_2122_2Ha_2Fa_ij =
    let s = Spin.Alfa in
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a ->
             h_two a k i n s s  *.  f_two a j l n s s 
          +. h_two a l i n s s  *.  f_two a j n k s s 
          -. h_two a k j n s s  *.  f_two a i l n s s   
          -. h_two a l j n s s  *.  f_two a i n k s s
        )
      )
  in

  if Parallel.master then Printf.printf "Computing m_2122_2Hb_2Fb_ij\n%!" ;
  let m_2122_2Hb_2Fb_ij =
    let s, s' = Spin.(Alfa, Beta) in
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a ->
            h_two a k n i s  s  *.  f_two a j n l s  s'
        +. h_two a l n j s  s' *.  f_two a i n k s  s 
        -. h_two a k j n s  s' *.  f_two a i l n s  s'
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_2122_2Hb_2Fb_ij2\n%!" ;
  let m_2122_2Hb_2Fb_ij2 =
    let s, s' = Spin.(Alfa, Beta) in
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
    sum mos_cabs (fun a ->
      -. h_two a l i n s  s' *.  f_two a j k n s  s' +.
      (if n < j then
                  h_two a k n i s  s' *.  f_two a j n l s' s'
               +. h_two a l n j s' s' *.  f_two a i n k s  s'
       else 
               -. h_two a k n i s  s' *.  f_two j a n l s' s' 
               -. h_two a l j n s' s' *.  f_two i a k n s  s'
      )
      )
    )
  in

  if Parallel.master then Printf.printf "Computing m_2122_2Ha_2Fa_ij2\n%!" ;
  let m_2122_2Ha_2Fa_ij2 =
    let s, s' = Spin.(Alfa, Beta) in
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a -> h_two a k n i s s' *.  f_two a j n l s s') +.
      sum mos_cabs (fun a -> h_two a l n j s s' *.  f_two a i n k s s') -.
      sum mos_cabs (fun a -> h_two a l n i s s' *.  f_two a j n k s s') -.
      sum mos_cabs (fun a -> h_two a k n j s s' *.  f_two a i n l s s')
      )
  in

  if Parallel.master then Printf.printf "Computing m_2122_2Ha_2Fa_nv\n%!" ;
  let m_2122_2Ha_2Fa_nv =
    let s = Spin.Alfa in
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a -> h_two a n i j s  s  *.  f_two a n l k s  s ) ) 
  in

  if Parallel.master then Printf.printf "Computing m_2122_2Hb_2Fb_nv\n%!" ;
  let m_2122_2Hb_2Fb_nv =
    let s, s' = Spin.(Alfa, Beta) in
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a -> -. h_two a n i j s s'  *.  f_two a n k l s s' ) ) 
  in

  if Parallel.master then Printf.printf "Computing m_2122_2Hb_2Fb_nv2\n%!" ;
  let m_2122_2Hb_2Fb_nv2 =
    let s, s' = Spin.(Alfa, Beta) in
    array_5_init mo_num mo_num mo_num mo_num mo_num (fun n i j k l ->
      sum mos_cabs (fun a -> -. h_two a n j i s s'  *.  f_two a n l k s s' ) ) 
  in

  let f_2 ki kj =
      let i, j, k, l, s, s', phase =
            match Excitation.of_det ki kj with
            | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
              hole, hole', particle, particle', spin, spin', phase
            | _ -> assert false
      in
      let mos_i, mos_i' =
        match s with
        | Spin.Alfa -> mos_a ki, mos_b ki
        | Spin.Beta -> mos_b ki, mos_a ki
      in

      let mos_j, mos_j' =
        match s with
        | Spin.Alfa -> mos_a kj, mos_b kj
        | Spin.Beta -> mos_b kj, mos_a kj
      in

      let mos_ij, mos_ij' =
        let alfa =
          let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
          let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
          Bitstring.to_list (Bitstring.logand i j)
          |> Array.of_list 
        in
        let beta  =
          let i = Spindeterminant.bitstring @@ Determinant.beta ki in
          let j = Spindeterminant.bitstring @@ Determinant.beta kj in
          Bitstring.to_list (Bitstring.logand i j)
          |> Array.of_list 
        in
        match s with
        | Spin.Alfa -> alfa, beta
        | Spin.Beta -> beta, alfa
      in

      let mos_novirt, mos_novirt' =
        let alfa =
          let i = Spindeterminant.bitstring @@ Determinant.alfa ki in
          let j = Spindeterminant.bitstring @@ Determinant.alfa kj in
          Bitstring.to_list (Bitstring.logor i j)
          |> Array.of_list 
        in
        let beta  =
          let i = Spindeterminant.bitstring @@ Determinant.beta ki in
          let j = Spindeterminant.bitstring @@ Determinant.beta kj in
          Bitstring.to_list (Bitstring.logor i j)
          |> Array.of_list 
        in
        match s with
        | Spin.Alfa -> alfa, beta
        | Spin.Beta -> beta, alfa
      in

      let result =
        if s = s' then
          m_2112_1H_2Fa.{i,j,k,l} +.
          sum mos_i      (fun n -> m_2112_2Ha_2Fa.{n,i,j,k,l} ) +.
          sum mos_i'     (fun n -> m_2112_2Hb_2Fa.{n,i,j,k,l} ) +.
          sum mos_j      (fun n -> m_2121_2Ha_2Fa.{n,i,j,k,l} ) +.
          sum mos_j'     (fun n -> m_2121_2Ha_2Fb.{n,i,j,k,l} ) +.
          sum mos_ij     (fun n -> m_2122_2Ha_2Fa_ij.{n,i,j,k,l} ) +.
          sum mos_ij'    (fun n -> m_2122_2Ha_2Fa_ij2.{n,i,j,k,l} ) +.
          sum mos_novirt (fun n -> m_2122_2Ha_2Fa_nv.{n,i,j,k,l} )
        else
          m_2112_1H_2Fb.{i,j,k,l} +.
          sum mos_i      (fun n -> m_2112_2Ha_2Fb.{n,i,j,k,l} ) +.
          sum mos_i'     (fun n -> m_2112_2Ha_2Fb.{n,j,i,l,k} ) +.
          sum mos_j      (fun n -> m_2121_2Hb_2Fa.{n,i,j,k,l} ) +.
          sum mos_j'     (fun n -> m_2121_2Hb_2Fa.{n,j,i,l,k} ) +.
          sum mos_novirt'(fun n -> m_2122_2Hb_2Fb_nv.{n,i,j,k,l} ) +.
          sum mos_novirt (fun n -> m_2122_2Hb_2Fb_nv2.{n,i,j,k,l} )+.
          sum mos_ij     (fun n -> m_2122_2Hb_2Fb_ij.{n,i,j,k,l} ) +.
          sum mos_ij'    (fun n -> m_2122_2Hb_2Fb_ij2.{n,i,j,k,l} )
      in

      match phase with
      | Phase.Pos -> result
      | Phase.Neg -> -. result

  in

  let ncabs = Array.length mos_cabs in
  let tmp_h = Mat.make0 ncabs 9 in
  let tmp_f = Mat.make0 ncabs 9 in

  let f_3 ki kj =
      let i, j, m, k, l, n, s1, s2, s3, phase =
            match Excitation.of_det ki kj with
            | Excitation.(Triple (phase,
                { hole=h1 ; particle=p1 ; spin=s1 },
                { hole=h2 ; particle=p2 ; spin=s2 },
                { hole=h3 ; particle=p3 ; spin=s3 }) ) ->
                      h1, h2, h3, p1, p2, p3, s1, s2, s3, phase
            | _ -> assert false
      in

      let result =
        let open Spin in
        match s1, s2, s3 with
        | Alfa, Alfa, Alfa
        | Beta, Beta, Beta ->
        (*
            sum mos_cabs (fun a -> h_two a k i j s1 s2 *. f_two a m n l s3 s3)
         +. sum mos_cabs (fun a -> h_two a n i j s1 s2 *. f_two a m l k s2 s3)
         -. sum mos_cabs (fun a -> h_two a l i j s1 s2 *. f_two a m n k s3 s3)
         +. sum mos_cabs (fun a -> h_two a l i m s1 s3 *. f_two a j n k s3 s2)
         -. sum mos_cabs (fun a -> h_two a n i m s1 s3 *. f_two a j l k s2 s2)
         -. sum mos_cabs (fun a -> h_two a k i m s1 s3 *. f_two a j n l s3 s2)
         +. sum mos_cabs (fun a -> h_two a n j m s2 s3 *. f_two a i l k s2 s1)
         +. sum mos_cabs (fun a -> h_two a k j m s2 s3 *. f_two a i n l s3 s1)
         -. sum mos_cabs (fun a -> h_two a l j m s2 s3 *. f_two a i n k s3 s1)

            sum_schwarz k i j m n l s1 s2 s3 s3
         +. sum_schwarz n i j m l k s1 s2 s2 s3
         -. sum_schwarz l i j m n k s1 s2 s3 s3
         +. sum_schwarz l i m j n k s1 s3 s3 s2
         -. sum_schwarz n i m j l k s1 s3 s2 s2
         -. sum_schwarz k i m j n l s1 s3 s3 s2
         +. sum_schwarz n j m i l k s2 s3 s2 s1
         +. sum_schwarz k j m i n l s2 s3 s3 s1
         -. sum_schwarz l j m i n k s2 s3 s3 s1
        *)
          Array.iteri (fun idx a -> tmp_f.{idx+1,7} <- f_two a i l k s2 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,5} <- f_two a j l k s2 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,2} <- f_two a m l k s2 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,9} <- f_two a i n k s3 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,4} <- f_two a j n k s3 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,3} <- f_two a m n k s3 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,8} <- f_two a i n l s3 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,6} <- f_two a j n l s3 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,1} <- f_two a m n l s3 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,3} <- if tmp_f.{idx+1,3} = 0. then 0. else  -. h_two a l i j s1 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,1} <- if tmp_f.{idx+1,1} = 0. then 0. else     h_two a k i j s1 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,2} <- if tmp_f.{idx+1,2} = 0. then 0. else     h_two a n i j s1 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,5} <- if tmp_f.{idx+1,5} = 0. then 0. else  -. h_two a n i m s1 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,6} <- if tmp_f.{idx+1,6} = 0. then 0. else  -. h_two a k i m s1 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,4} <- if tmp_f.{idx+1,4} = 0. then 0. else     h_two a l i m s1 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,9} <- if tmp_f.{idx+1,9} = 0. then 0. else  -. h_two a l j m s2 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,7} <- if tmp_f.{idx+1,7} = 0. then 0. else     h_two a n j m s2 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,8} <- if tmp_f.{idx+1,8} = 0. then 0. else     h_two a k j m s2 s3) mos_cabs;
          dot (Mat.as_vec tmp_h) (Mat.as_vec tmp_f)
        | Alfa, Alfa, Beta
        | Beta, Beta, Alfa ->
        (*
            sum mos_cabs (fun a -> h_two a l i j s1 s2 *. f_two a m k n s1 s3)
         -. sum mos_cabs (fun a -> h_two a k i j s1 s2 *. f_two a m l n s1 s3)
         +. sum mos_cabs (fun a -> h_two a l m j s3 s2 *. f_two a i n k s3 s1)
         -. sum mos_cabs (fun a -> h_two a k m j s3 s2 *. f_two a i n l s3 s1)
         +. sum mos_cabs (fun a -> h_two a k m i s3 s1 *. f_two a j n l s3 s2)
         -. sum mos_cabs (fun a -> h_two a l m i s3 s1 *. f_two a j n k s3 s2)
         +. sum mos_cabs (fun a -> h_two a n j m s2 s3 *. f_two a i l k s2 s1)
         -. sum mos_cabs (fun a -> h_two a n i m s1 s3 *. f_two a j l k s2 s2)

            sum_schwarz l i j m k n s1 s2 s1 s3
         -. sum_schwarz k i j m l n s1 s2 s1 s3
         +. sum_schwarz l m j i n k s3 s2 s3 s1
         -. sum_schwarz k m j i n l s3 s2 s3 s1
         +. sum_schwarz k m i j n l s3 s1 s3 s2
         -. sum_schwarz l m i j n k s3 s1 s3 s2
         +. sum_schwarz n j m i l k s2 s3 s2 s1
         -. sum_schwarz n i m j l k s1 s3 s2 s2
        *)
          Array.iteri (fun idx a -> tmp_f.{idx+1,1} <- f_two a m k n s1 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,2} <- f_two a m l n s1 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,4} <- f_two a i n l s3 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,5} <- f_two a j n l s3 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,3} <- f_two a i n k s3 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,6} <- f_two a j n k s3 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,7} <- f_two a i l k s2 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,8} <- f_two a j l k s2 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,1} <- if tmp_f.{idx+1,1} = 0. then 0. else    h_two a l i j s1 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,2} <- if tmp_f.{idx+1,2} = 0. then 0. else -. h_two a k i j s1 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,3} <- if tmp_f.{idx+1,3} = 0. then 0. else +. h_two a l m j s3 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,4} <- if tmp_f.{idx+1,4} = 0. then 0. else -. h_two a k m j s3 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,5} <- if tmp_f.{idx+1,5} = 0. then 0. else +. h_two a k m i s3 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,6} <- if tmp_f.{idx+1,6} = 0. then 0. else -. h_two a l m i s3 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,7} <- if tmp_f.{idx+1,7} = 0. then 0. else +. h_two a n j m s2 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,8} <- if tmp_f.{idx+1,8} = 0. then 0. else -. h_two a n i m s1 s3) mos_cabs;
          dot ~n:(8*ncabs) (Mat.as_vec tmp_h) (Mat.as_vec tmp_f)
        | Alfa, Beta, Beta
        | Beta, Alfa, Alfa ->
        (*
            sum mos_cabs (fun a -> h_two a l i j s1 s2 *. f_two a m k n s1 s3)
         -. sum mos_cabs (fun a -> h_two a n i j s1 s2 *. f_two a m k l s1 s2)
         +. sum mos_cabs (fun a -> h_two a n i m s1 s3 *. f_two a j k l s1 s2)
         +. sum mos_cabs (fun a -> h_two a n j m s2 s3 *. f_two a i l k s2 s1)
         -. sum mos_cabs (fun a -> h_two a l i m s1 s3 *. f_two a j k n s1 s3)
         -. sum mos_cabs (fun a -> h_two a l j m s2 s3 *. f_two a i n k s3 s1)
         +. sum mos_cabs (fun a -> h_two a k m i s3 s1 *. f_two a j n l s3 s2)
         -. sum mos_cabs (fun a -> h_two a k j i s2 s1 *. f_two a m n l s3 s2)

            sum_schwarz l i j m k n s1 s2 s1 s3
         -. sum_schwarz n i j m k l s1 s2 s1 s2
         +. sum_schwarz n i m j k l s1 s3 s1 s2
         +. sum_schwarz n j m i l k s2 s3 s2 s1
         -. sum_schwarz l i m j k n s1 s3 s1 s3
         -. sum_schwarz l j m i n k s2 s3 s3 s1
         +. sum_schwarz k m i j n l s3 s1 s3 s2
         -. sum_schwarz k j i m n l s2 s1 s3 s2
         *)
          Array.iteri (fun idx a -> tmp_f.{idx+1,6} <-    f_two a i n k s3 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,4} <-    f_two a i l k s2 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,5} <-    f_two a j k n s1 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,1} <-    f_two a m k n s1 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,2} <-    f_two a m k l s1 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,3} <-    f_two a j k l s1 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,7} <-    f_two a j n l s3 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_f.{idx+1,8} <-    f_two a m n l s3 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,1} <- if tmp_f.{idx+1,1} = 0. then 0. else    h_two a l i j s1 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,2} <- if tmp_f.{idx+1,2} = 0. then 0. else -. h_two a n i j s1 s2) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,3} <- if tmp_f.{idx+1,3} = 0. then 0. else +. h_two a n i m s1 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,4} <- if tmp_f.{idx+1,4} = 0. then 0. else +. h_two a n j m s2 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,5} <- if tmp_f.{idx+1,5} = 0. then 0. else -. h_two a l i m s1 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,6} <- if tmp_f.{idx+1,6} = 0. then 0. else -. h_two a l j m s2 s3) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,7} <- if tmp_f.{idx+1,7} = 0. then 0. else +. h_two a k m i s3 s1) mos_cabs;
          Array.iteri (fun idx a -> tmp_h.{idx+1,8} <- if tmp_f.{idx+1,8} = 0. then 0. else -. h_two a k j i s2 s1) mos_cabs;
          dot ~n:(8*ncabs) (Mat.as_vec tmp_h) (Mat.as_vec tmp_f)
        | Beta, Alfa, Beta
        | Alfa, Beta, Alfa ->  assert false
      in
      match phase with
      | Phase.Pos -> result
      | Phase.Neg -> -. result
  in

  let result = 
    { simulation ; aux_basis ;
      f_0 ; f_1 ; f_2 ; f_3 
    } 
  in
  result