diff --git a/run_test_f12.ml b/run_test_f12.ml
new file mode 100644
index 0000000..d73d625
--- /dev/null
+++ b/run_test_f12.ml
@@ -0,0 +1,2309 @@
+
+let png_image = print_endline ;;
+
+
+open Lacaml.D
+
+
+let basis_filename  = "/home/scemama/qp2/data/basis/6-31g" 
+let aux_basis_filename = "/home/scemama/qp2/data/basis/cc-pvdz" 
+let nuclei = Nuclei.of_zmt_string "li" 
+let frozen_core = false
+let multiplicity = 2
+let state = 1
+
+let f12 = F12factor.gaussian_geminal 1.0 
+let basis =  Basis.of_nuclei_and_basis_filenames ~f12 ~nuclei [basis_filename] 
+let charge = 0 
+
+
+let simulation =
+    Simulation.make 
+        ~charge ~multiplicity ~nuclei
+        ~cartesian:true
+        basis
+
+
+let n_elec_alfa, n_elec_beta, n_elec = 
+  let e =  Simulation.electrons simulation in
+  Electrons.(n_alfa e, n_beta e, n_elec e)
+
+let hf =  HartreeFock.make ~guess:`Hcore ~max_scf:1 simulation  ;;
+
+let mo_basis = MOBasis.of_hartree_fock hf
+
+let f12 = Util.of_some @@ Simulation.f12 simulation 
+
+let mo_num = MOBasis.size mo_basis 
+
+let pp_spindet = Spindeterminant.pp mo_num
+
+let pp_det = Determinant.pp mo_num
+
+;;
+
+
+
+let simulation_aux = 
+    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 
+
+
+let aux_basis = 
+  MOBasis.of_mo_basis simulation_aux mo_basis
+
+let aux_num = 
+  MOBasis.size aux_basis
+
+
+let () = ignore @@ MOBasis.f12_ints   aux_basis
+let () = ignore @@ MOBasis.two_e_ints aux_basis
+
+let cancel_singles = false 
+
+let mos_cabs = 
+  Util.list_range (mo_num+1) aux_num
+  
+let mos_abs = 
+  Util.list_range 1 aux_num
+  
+let mos_in = 
+  Util.list_range 1 mo_num
+
+let mos_a k =
+  Determinant.alfa k
+  |> Spindeterminant.to_list
+  
+let mos_b k =
+  Determinant.beta k
+  |> Spindeterminant.to_list
+  
+
+let h_one =
+  let h = 
+    MOBasis.one_e_ints aux_basis
+  in fun i j _ -> h.{i,j}
+  
+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
+
+    
+  
+
+let f_two = 
+  let two_e_ints = MOBasis.f12_ints aux_basis in
+  let f2 i j k l (s:Spin.t) (s':Spin.t) =
+    if s' <> s then
+        0.5 *. F12.get_phys two_e_ints i j k l
+    else
+        0.5 *. (
+         (F12.get_phys two_e_ints i j k l) -.                                                 
+         (F12.get_phys two_e_ints i j l k) )
+  in
+  let f3 i j k l (s:Spin.t) (s':Spin.t) = 
+      if (i=k && j<>l) || (j=l && i<>k) then
+        0.
+      else
+        f2 i j k l s s'
+  in
+  if cancel_singles then f3 else f2
+
+let f_one = fun _ _ _ -> 0.
+
+(*
+let f_two = h_two
+
+let f_one = h_one
+*)
+
+let f12_integrals mo_basis =
+  ( f_one, f_two, None )
+
+let h_ij mo_basis ki kj =
+  let integrals =
+    List.map (fun f -> f mo_basis)
+      [ CI.h_integrals ]
+  in
+  CIMatrixElement.make integrals ki kj 
+  |> List.hd
+
+
+let f_ij mo_basis ki kj =
+  let integrals =
+    List.map (fun f -> f mo_basis)
+      [ f12_integrals ]
+  in
+  CIMatrixElement.make integrals ki kj 
+  |> List.hd                                                                                  
+
+
+let hf_ij mo_basis ki kj =
+  let integrals =
+    List.map (fun f -> f mo_basis)
+      [ CI.h_integrals ; f12_integrals ]
+  in
+  CIMatrixElement.make integrals ki kj
+
+
+
+
+let is_a_double det_space =
+    let mo_class = DeterminantSpace.mo_class det_space in
+    let mo_num = Array.length @@ MOClass.mo_class_array mo_class in
+    let m l =
+      List.fold_left (fun accu i ->
+          let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
+        ) (Bitstring.zero mo_num) l
+    in
+    let aux_mask   = m (MOClass.auxiliary_mos mo_class) in
+    fun k ->                                                                                  
+      let alfa =
+        Determinant.alfa k
+        |> Spindeterminant.bitstring
+      in
+      let beta =
+        Determinant.beta k
+        |> Spindeterminant.bitstring
+      in
+      let a = Bitstring.logand aux_mask alfa
+      and b = Bitstring.logand aux_mask beta
+      in
+      match Bitstring.popcount a + Bitstring.popcount b with
+      | 2 | 1 -> true
+      | 0 | _ -> false
+
+
+
+let in_space = 
+    DeterminantSpace.fci_f12_of_mo_basis aux_basis ~frozen_core mo_num
+    
+let aux_space = 
+    DeterminantSpace.fci_of_mo_basis aux_basis ~frozen_core 
+
+let det_space_in () =
+    DeterminantSpace.determinant_stream in_space
+
+let det_space_out () =
+    let s = 
+        DeterminantSpace.determinant_stream aux_space
+    in
+    Stream.from (fun _ ->
+        try
+          let is_a_double = is_a_double in_space in
+          let rec result () =
+            let ki = Stream.next s in
+            if is_a_double ki then
+              Some (ki,ki)
+            else
+              result ()
+          in
+          result ()
+        with Stream.Failure -> None
+      )
+    
+
+let ci = CI.make ~n_states:state in_space
+
+let ci_coef, ci_energy = Lazy.force ci.eigensystem 
+
+let _ = print_newline () 
+
+let p12 det_space =
+    let mo_class = DeterminantSpace.mo_class det_space in
+    let mo_num = Array.length @@ MOClass.mo_class_array mo_class in
+    let m l =
+      List.fold_left (fun accu i ->
+          let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
+        ) (Bitstring.zero mo_num) l
+    in
+    let aux_mask = m (MOClass.auxiliary_mos mo_class) in
+    let not_aux_mask =
+      Bitstring.(shift_left_one mo_num (mo_num-1) |> minus_one |> logxor aux_mask)
+    in
+    fun k ->
+      let alfa =
+        Determinant.alfa k
+        |> Spindeterminant.bitstring
+      in
+      let beta =
+        Determinant.beta k
+        |> Spindeterminant.bitstring
+      in
+      let a = Bitstring.logand aux_mask alfa
+      and b = Bitstring.logand aux_mask beta
+      in
+      match Bitstring.popcount a, Bitstring.popcount b with
+      | 2, 0  
+      | 0, 2 -> Some (Determinant.negate_phase k)
+      | 1, 1 -> Some (Determinant.of_spindeterminants
+                  (Spindeterminant.of_bitstring  @@
+                    Bitstring.(logor b (logand not_aux_mask alfa)) )
+                  (Spindeterminant.of_bitstring  @@
+                    Bitstring.(logor a (logand not_aux_mask beta))
+                  ) )
+      | 1, 0 
+      | 0, 1 -> Some k
+      | _ -> None
+                    
+        
+
+let out_list =
+    Util.stream_to_list (det_space_out ())
+    
+let in_list =
+    Util.stream_to_list (det_space_in ())
+    
+let det_a = Array.of_list out_list
+            |> Array.map (fun (i,_) -> i)
+
+let det_I =
+    let n = 123456789 in
+    in_list
+    |> List.map (fun k -> (Random.int n, k))
+    |> List.sort compare
+    |> List.map  (fun (_,k) -> k)
+    |> Array.of_list
+
+
+let generate_alphas ki kj exc cabs l r =
+  (* Check input excitation degree *)
+  let d = Determinant.degree ki kj in
+  if d <> exc then
+    Printf.sprintf "Invalid excitation degree. Expected %d and found %d." exc d
+    |> failwith;
+    
+  (* Generate single excitations *)
+  let all_singles ki =
+      let mos_a, mos_b = Determinant.to_lists ki in
+      [ List.map (fun hole -> 
+          List.map (fun particle ->
+               if hole = particle then None else
+                  Some (Determinant.single_excitation Spin.Alfa hole particle ki)
+              ) mos_abs
+          ) mos_a
+      ;
+        List.map (fun hole -> 
+          List.map (fun particle ->
+               if hole = particle then None else
+                 Some (Determinant.single_excitation Spin.Beta hole particle ki)
+              ) mos_abs
+          ) mos_b 
+      ]
+      |> List.concat
+      |> List.concat
+      |> Util.list_some
+      |> List.filter (fun x -> Determinant.is_none x = false)
+      |> List.map (fun x -> Determinant.set_phase Phase.Pos x)
+  in
+  
+  (* Generate double excitations *)
+  let all_doubles ki =
+      let mos_a, mos_b = Determinant.to_lists ki in
+      [ List.map (fun hole ->  (* Alpha-Alpha *)
+          List.map (fun particle ->
+            List.map (fun hole' -> 
+              List.map (fun particle' ->
+                if hole = particle || hole' = particle' then None else
+                if hole' > hole && particle' < particle then
+                  Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Alfa hole' particle' ki)
+                else None
+                      ) mos_abs
+                  ) mos_a 
+              ) mos_abs
+          ) mos_a
+        ; 
+        List.map (fun hole ->  (* Beta-Beta *)
+          List.map (fun particle ->
+            List.map (fun hole' -> 
+              List.map (fun particle' ->
+                if hole' > hole && particle' < particle && hole <> particle && hole' <> particle' then
+                  Some (Determinant.double_excitation Spin.Beta hole particle Spin.Beta hole' particle' ki)
+                else None
+                      ) mos_abs
+                  ) mos_b 
+              ) mos_abs
+          ) mos_b
+        ;
+        List.map (fun hole ->  (* Alpha-Beta *)
+          List.map (fun particle ->
+            List.map (fun hole' -> 
+              List.map (fun particle' ->
+                  if hole = particle || hole' = particle' then None else
+                  Some (Determinant.double_excitation Spin.Alfa hole particle Spin.Beta hole' particle' ki)
+                      ) mos_abs
+                  ) mos_b 
+              ) mos_abs
+          ) mos_a
+      ]
+      |> List.concat
+      |> List.concat
+      |> List.concat
+      |> List.concat
+      |> Util.list_some
+      |> List.filter (fun x -> Determinant.is_none x = false)
+      |> List.map (fun x -> Determinant.set_phase Phase.Pos x)
+  in
+
+  (* Generate left and right excitations *)
+  let al = 
+      match l with
+      | 1 -> all_singles ki
+      | 2 -> all_doubles ki
+      | _ -> assert false
+  in
+  let ar = 
+      match r with
+      | 1 -> all_singles kj
+      | 2 -> all_doubles kj
+      | _ -> assert false
+  in
+  
+  let mo_class = DeterminantSpace.mo_class in_space in
+  let m l =
+    List.fold_left (fun accu i ->
+        let j = i-1 in Bitstring.logor accu (Bitstring.shift_left_one mo_num j)
+      ) (Bitstring.zero mo_num) l
+  in
+  let aux_mask = m (MOClass.auxiliary_mos mo_class) in
+  let good_cabs k = 
+      let alfa =
+        Determinant.alfa k
+        |> Spindeterminant.bitstring
+      in
+      let beta =
+        Determinant.beta k
+        |> Spindeterminant.bitstring
+      in
+      let a = Bitstring.logand aux_mask alfa
+      and b = Bitstring.logand aux_mask beta
+      in
+      Bitstring.popcount a + Bitstring.popcount b = cabs
+  in
+  let good_lr k =
+    Determinant.degree ki k = l &&
+    Determinant.degree k kj = r 
+  in
+
+
+  
+  (* Merge lists in a set of unique determinants *)
+  List.concat [ al; ar ]
+  |> List.sort_uniq compare
+  
+  (* Filter out all determinants with incorrect numbers of electrons in the CABS *)
+  |> List.filter good_cabs 
+  
+  (* Filter out all determinants with incorrect excitation numbers *)
+  |> List.filter good_lr 
+  
+  
+let compute_HaaF ki alphas kj =
+   List.fold_left (fun accu alpha -> accu
+      +. h_ij aux_basis ki alpha
+      *. f_ij aux_basis alpha kj
+    ) 0. alphas
+
+
+let check n integral_value exc cabs lexc rexc =
+  let cpudet, cpuint = ref 0., ref 0. in
+  let det_list =
+      match n with
+        | 0 -> det_I
+        | n -> Array.sub det_I 0 n
+  in
+  let result =
+        if Parallel.master then Printf.printf "Checking ... \n%!";
+        let percent = ref 0 in
+        let task (i,ki) =
+            (let p = (10 * (i+1))/(Array.length det_list) in
+            if p <> !percent then
+                ( percent := p ; if Parallel.master then Printf.printf " - %3d %%\n%!" (p*10) ));
+            Array.mapi (fun j kj -> 
+                if i > j || Determinant.degree ki kj <> exc then
+                    (0,0,0.,0.)
+                else
+                  begin
+                    let alphas = generate_alphas ki kj exc cabs lexc rexc in
+                    let det_value =
+                      let t0 = Unix.gettimeofday () in
+                      let result = compute_HaaF ki alphas kj in
+                      cpudet := !cpudet +. Unix.gettimeofday () -. t0;
+                      result
+                    in
+                    let int_value = 
+                      let t0 = Unix.gettimeofday () in
+                      let result = integral_value ki kj in
+                      cpuint := !cpuint +. Unix.gettimeofday () -. t0;
+                      result
+                    in
+(*   Printf.printf "%d %d %e %e\n%!" i j det_value int_value;  *)
+                    (i,j,det_value,int_value)
+                  end
+            ) det_list
+        in
+        det_list
+        |> Array.mapi (fun i ki -> (i,ki))
+        |> Array.to_list
+        |> Stream.of_list
+        |> Farm.run ~f:task
+        |> Util.stream_to_list
+        |> Array.concat
+  in
+  let i,j,d,v = 
+      let rec aux k imax jmax emax dmax vmax = 
+        if k = -1 then
+          imax, jmax, dmax, vmax
+        else
+          let i, j, d, v = result.(k) in
+          let e = abs_float (d -. v) in
+          if e >= emax then
+            aux (k-1) i j e d v
+          else
+            aux (k-1) imax jmax emax dmax vmax
+      in
+      aux (Array.length result - 1) 0 0 0. 0. 0.
+  in
+  let error = abs_float (d -. v) in
+  if error < epsilon_float then
+  (*
+      Printf.printf "OK: %e\n%!" error
+      *)
+      if Parallel.master then Printf.printf "OK: (%d, %d) | %e %e | %e  |  cpu : %f  %f\n%!" i j d v error !cpudet !cpuint
+  else
+      if Parallel.master then Printf.printf "Failed: (%d, %d) | %e %e | %e  |  cpu : %f %f\n%!" i j d v error !cpudet !cpuint
+
+    
+let sum l f = List.fold_left (fun accu i -> accu +. f i) 0. l
+
+let array_3_init d1 d2 d3 f =
+  let result = 
+      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
+             result.{i,j,k} <- f i j k
+          done
+      done
+  done;
+  result
+
+
+let array_4_init d1 d2 d3 d4 f =
+  let result = 
+      Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4 |]
+  in
+  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} <- f i j k l
+        done
+      done
+    done
+  done;
+  result
+
+
+let array_5_init d1 d2 d3 d4 d5 f =
+  let result = 
+      Bigarray.(Genarray.create Float64 fortran_layout) [| d1;d2;d3;d4;d5 |]
+  in
+  for m=1 to d5 do
+    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} <- f i j k l m
+          done
+        done
+      done
+    done
+  done;
+  result
+
+(* ----- *)
+
+let _ = png_image "0_1_11.png"
+
+let integral_value ki kj = 
+   (* Alpha-Beta *)
+   let s, s' = Spin.(Alfa, Beta) in 
+   let mos_a, mos_b = mos_a ki, mos_b ki in
+   
+   sum mos_cabs (fun a ->
+      sum mos_a (fun i -> 
+           (h_one i a s  +.
+               sum mos_a (fun j -> h_two i j a j s s ) +.
+               sum mos_b (fun j -> h_two i j a j s s') ) *. 
+           (f_one a i s  +.
+               sum mos_a (fun j -> f_two a j i j s s ) +.
+               sum mos_b (fun j -> f_two a j i j s s') )
+         )
+       +.
+       sum mos_b (fun i -> 
+           (h_one i a s  +.
+               sum mos_b (fun j -> h_two i j a j s' s') +.
+               sum mos_a (fun j -> h_two i j a j s' s ) ) *. 
+           (f_one a i s  +.
+               sum mos_a (fun j -> f_two a j i j s' s ) +.
+               sum mos_b (fun j -> f_two a j i j s' s') )
+         ) 
+     ) 
+     
+let _ =
+  check 0 integral_value 0 1 1 1
+    
+
+let m_0111_1H_1F = 
+  Vec.init mo_num (fun i -> 
+        sum mos_cabs (fun a ->
+            h_one i a Spin.Alfa *. f_one a i Spin.Alfa ))
+
+    
+let m_0111_1H_2Fa = 
+  Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->
+    h_one i a Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Alfa  ))
+
+    
+let m_0111_1H_2Fb = 
+  Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->
+    h_one i a Spin.Alfa *. f_two a j i j Spin.Alfa Spin.Beta  ))
+
+
+let m_0111_2Ha_1F = 
+  Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->
+    h_two i j a j Spin.Alfa Spin.Alfa *. f_one a j Spin.Alfa ))
+
+    
+let m_0111_2Hb_1F = 
+  Mat.init_cols mo_num mo_num (fun i j -> sum mos_cabs (fun a ->
+    h_two i j a j Spin.Alfa Spin.Beta *. f_one a j Spin.Alfa ))
+
+
+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 i j a j Spin.Alfa Spin.Alfa *.
+      f_two a k i k Spin.Alfa Spin.Alfa 
+    )
+  )
+
+
+let m_0111_2Hb_2Fa =
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+    sum mos_cabs (fun a ->
+      h_two i j a j Spin.Alfa Spin.Beta *.
+      f_two a k i k Spin.Alfa Spin.Alfa 
+    )
+  )
+
+
+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 i j a j Spin.Alfa Spin.Alfa *.
+      f_two a k i k Spin.Alfa Spin.Beta 
+    )
+  )
+
+
+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 i j a j Spin.Alfa Spin.Beta *.
+      f_two a k i k Spin.Alfa Spin.Beta 
+    )
+  )
+
+
+
+let integral_value ki kj = 
+
+   let mos_a, mos_b = mos_a ki, mos_b ki in
+   
+   sum mos_a (fun i ->
+     m_0111_1H_1F.{i} +. 
+     sum mos_a (fun j ->
+       m_0111_1H_2Fa.{i,j} +. m_0111_2Ha_1F.{i,j} +.
+       sum mos_a (fun k -> m_0111_2Ha_2Fa.{i,j,k}) +.
+       sum mos_b (fun k -> m_0111_2Ha_2Fb.{i,j,k})) +.
+     sum mos_b (fun j ->
+       m_0111_1H_2Fb.{i,j} +. m_0111_2Hb_1F.{i,j} +.
+       sum mos_a (fun k -> m_0111_2Hb_2Fa.{i,j,k}) +.
+       sum mos_b (fun k -> m_0111_2Hb_2Fb.{i,j,k}))
+     )
+   +.
+   sum mos_b (fun i ->
+     m_0111_1H_1F.{i} +. 
+     sum mos_b (fun j ->
+       m_0111_1H_2Fa.{i,j} +. m_0111_2Ha_1F.{i,j} +.
+       sum mos_b (fun k -> m_0111_2Ha_2Fa.{i,j,k}) +.
+       sum mos_a (fun k -> m_0111_2Ha_2Fb.{i,j,k})) +.
+     sum mos_a (fun j ->
+       m_0111_1H_2Fb.{i,j} +. m_0111_2Hb_1F.{i,j} +.
+       sum mos_b (fun k -> m_0111_2Hb_2Fa.{i,j,k}) +.
+       sum mos_a (fun k -> m_0111_2Hb_2Fb.{i,j,k}))
+     )
+   
+let _ =
+  check 0 integral_value 0 1 1 1
+    
+
+let _ = png_image "0_1_22.png"
+
+let integral_value ki kj = 
+   (* mos unoccupied in both I and J *)
+   let mos_virt_a, mos_virt_b = 
+     Array.init mo_num (fun i -> Some (i+1)) , 
+     Array.init mo_num (fun i -> Some (i+1)) 
+   in
+   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
+   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
+   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
+   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
+   
+   let mos_virt_a, mos_virt_b = 
+     Array.to_list mos_virt_a |> Util.list_some,
+     Array.to_list mos_virt_b |> Util.list_some
+   in
+   
+   (* Alpha-Beta *)
+   let s, s' = Spin.(Alfa, Beta) in
+   sum mos_cabs   (fun a ->
+   sum (mos_b ki) (fun j ->  
+   sum (mos_a ki) (fun i ->  
+   sum mos_virt_b (fun k ->  
+       h_two i j a k s s'  *. f_two a k i j s s'
+   )))
+   +.
+   sum (mos_a ki) (fun j ->  
+   sum (mos_b ki) (fun i ->  
+   sum mos_virt_a (fun k ->  
+       h_two i j a k s s'  *. f_two a k i j s s'
+   )))) 
+   +.
+   (* Alpha-Alpha / Beta-Beta *)
+   List.fold_left (fun accu (mos_virt,mos,s) -> 
+     let s' = s in accu +.
+     sum mos (fun j ->  
+     sum mos (fun i -> if i < j then 0. else
+     sum mos_virt (fun k ->  
+     sum mos_cabs (fun a -> 
+       h_two i j a k s s'  *. f_two a k i j s s'
+     )))) 
+   ) 0. [ (mos_virt_a,mos_a ki,Spin.Alfa) ; (mos_virt_b,mos_b ki,Spin.Beta)]
+     
+let _ =
+  check 0 integral_value 0 1 2 2
+    
+
+let m_0122_Haa = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+      sum mos_cabs (fun a ->
+        h_two i j a k Spin.Alfa Spin.Alfa *. f_two a k i j Spin.Alfa Spin.Alfa
+        ) )
+      
+let m_0122_Hab = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+      sum mos_cabs (fun a ->
+         h_two i j a k Spin.Alfa Spin.Beta *. f_two a k i j Spin.Alfa Spin.Beta
+      ) )
+
+
+let integral_value ki kj = 
+   (* mos unoccupied in both I and J *)
+   let mos_virt_a, mos_virt_b = 
+     Array.init mo_num (fun i -> Some (i+1)) , 
+     Array.init mo_num (fun i -> Some (i+1)) 
+   in
+   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
+   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
+   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
+   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
+   
+   let mos_virt_a, mos_virt_b = 
+     Array.to_list mos_virt_a |> Util.list_some,
+     Array.to_list mos_virt_b |> Util.list_some
+   in
+   
+   let result = 
+   sum (mos_a ki) (fun i ->  
+   sum (mos_b ki) (fun j ->  
+   sum mos_virt_b (fun k ->  m_0122_Hab.{i,j,k} )))
+   +.
+   sum (mos_b ki) (fun i ->  
+   sum (mos_a ki) (fun j ->  
+   sum mos_virt_a (fun k ->  m_0122_Hab.{i,j,k} )))
+   +.
+   (* Alpha-Alpha / Beta-Beta *)
+   List.fold_left (fun accu (mos_virt,mos) -> 
+     accu +.
+     sum mos (fun j ->  
+     sum mos (fun i -> if i <= j then 0. else
+     sum mos_virt (fun k -> m_0122_Haa.{i,j,k} )))
+   ) 0. [ (mos_virt_a,mos_a ki) ; (mos_virt_b,mos_b ki)]
+   in 
+   result
+     
+let _ =
+  check 0 integral_value 0 1 2 2
+    
+
+let _ = png_image "0_2_22.png"
+
+let integral_value ki kj = 
+   (* Alpha-Beta *)
+   let s, s' = Spin.(Alfa, Beta) in
+   sum  mos_cabs  (fun a -> 
+   sum  mos_cabs  (fun b ->  
+   sum (mos_a ki) (fun i ->  
+   sum (mos_b ki) (fun j ->  
+       h_two i j a b s s'  *. f_two a b i j s s'
+   ))))
+   +.
+   (* Alpha-Alpha / Beta-Beta *)
+   List.fold_left (fun accu (mos,s) -> 
+     let s' = s in accu +.
+     sum mos_cabs (fun b ->  
+     sum mos_cabs (fun a -> if b > a then 0. else
+     sum   mos    (fun j ->  
+     sum   mos    (fun i -> if i < j then 0. else
+        h_two i j a b s s'  *. f_two a b i j s s'
+     ))))
+   ) 0. [ (mos_a ki,Spin.Alfa) ; (mos_b ki,Spin.Beta)]
+     
+let _ =
+  check 0 integral_value 0 2 2 2
+    
+
+let m_0222_Haa = 
+  Mat.init_cols mo_num mo_num (fun i j ->
+      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
+        ) )
+     )
+        
+let m_0222_Hab = 
+  Mat.init_cols mo_num mo_num (fun i j ->
+      sum mos_cabs (fun a ->
+      sum mos_cabs (fun b ->
+        h_two i j a b Spin.Alfa Spin.Beta *. f_two a b i j Spin.Alfa Spin.Beta
+        ) )
+     )
+      
+let integral_value ki kj = 
+   (* Alpha-Beta *)
+   sum (mos_a ki) (fun i ->  
+       sum (mos_b ki) (fun j ->  
+          m_0222_Hab.{i,j}  
+       ))
+   +.
+   (* Alpha-Alpha / Beta-Beta *)
+   List.fold_left (fun accu (mos) -> 
+     accu +.
+     sum mos (fun j ->  
+         sum mos (fun i -> if i < j then 0. else
+           m_0222_Haa.{i,j}  
+     ))
+   ) 0. [ (mos_a ki) ; (mos_b ki)]
+     
+let _ =
+  check 0 integral_value 0 2 2 2
+    
+
+let _ = png_image "1_1_11.png"
+
+let integral_value 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_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 = 
+       let s' = Spin.other s in
+       
+       sum mos_cabs (fun a -> 
+           (h_one i a s +. 
+              sum mos_i  (fun j ->  h_two i j a j s s ) +.
+              sum mos_i' (fun j ->  h_two i j a j s s') ) *.
+           (f_one a k s +. 
+              sum mos_j  (fun j ->  f_two a j k j s s ) +.
+              sum mos_j' (fun j ->  f_two a j k j s s') )  
+         )
+    in
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 1 1 1 1
+
+let m_1111_1H_1F = 
+  Mat.init_cols mo_num mo_num (fun i k -> 
+        sum mos_cabs (fun a ->
+            h_one i a Spin.Alfa *. f_one a k Spin.Alfa ))
+
+
+let m_1111_1H_2Fa = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+      sum mos_cabs (fun a ->
+        h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Alfa ))
+
+    
+let m_1111_1H_2Fb = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+      sum mos_cabs (fun a ->
+        h_one i a Spin.Alfa *. f_two a j k j Spin.Alfa Spin.Beta ))
+
+
+let m_1111_2Ha_1F = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+   sum mos_cabs (fun a ->
+    h_two i j a j Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa ))
+
+    
+let m_1111_2Hb_1F = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+   sum mos_cabs (fun a ->
+    h_two i j a j Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa ))
+
+
+let m_1111_2Ha_2Fa =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+    sum mos_cabs (fun a ->
+      h_two i j a j Spin.Alfa Spin.Alfa *.
+      f_two a l k l Spin.Alfa Spin.Alfa 
+    )
+  )
+
+
+let m_1111_2Hb_2Fa =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+    sum mos_cabs (fun a ->
+      h_two i j a j Spin.Alfa Spin.Beta *.
+      f_two a l k l Spin.Alfa Spin.Alfa 
+    )
+  )
+
+
+let m_1111_2Ha_2Fb =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+    sum mos_cabs (fun a ->
+      h_two i j a j Spin.Alfa Spin.Alfa *.
+      f_two a l k l Spin.Alfa Spin.Beta 
+    )
+  )
+
+
+let m_1111_2Hb_2Fb =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+    sum mos_cabs (fun a ->
+      h_two i j a j Spin.Alfa Spin.Beta *.
+      f_two a l k l Spin.Alfa Spin.Beta 
+    )
+  )
+
+
+let integral_value 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_a, mos_b =
+      match s with
+      | Spin.Alfa -> mos_a ki, mos_b ki
+      | Spin.Beta -> mos_b ki, mos_a ki
+    in    
+   
+    let result = 
+     m_1111_1H_1F.{i,k} +. 
+     sum mos_a (fun j -> if j = i then 0. else
+       m_1111_1H_2Fa.{i,j,k} +. m_1111_2Ha_1F.{i,j,k} +.
+       sum mos_a (fun l -> if l = i then 0. else m_1111_2Ha_2Fa.{i,j,k,l}) +.
+       sum mos_b (fun l -> m_1111_2Ha_2Fb.{i,j,k,l})) +.
+     sum mos_b (fun j ->
+       m_1111_1H_2Fb.{i,j,k} +. m_1111_2Hb_1F.{i,j,k} +.
+       sum mos_a (fun l -> if l = i then 0. else m_1111_2Hb_2Fa.{i,j,k,l}) +.
+       sum mos_b (fun l -> m_1111_2Hb_2Fb.{i,j,k,l}))
+    in
+    
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+   
+let _ =
+  check 0 integral_value 1 1 1 1
+    
+
+let _ = png_image "1_1_12.png"
+
+let integral_value 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_j, mos_j' =
+      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)
+      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)
+      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 result = 
+       let s' = Spin.other s in
+       sum mos_cabs (fun a -> 
+           sum mos_j (fun j -> 
+               (h_one j a s +. 
+                  sum mos_i  (fun m -> h_two j m a m s s )  +.
+                  sum mos_i' (fun m -> h_two j m a m s s')  ) *.
+               (f_two a i j k s s )
+             ) 
+           +.
+           sum mos_j' (fun j -> 
+               (h_one j a s +. 
+                  sum mos_i  (fun m ->  h_two j m a m s' s ) +.
+                  sum mos_i' (fun m ->  h_two j m a m s' s') ) *.
+               (f_two a i j k s' s )
+             ) 
+         ) 
+    in
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 1 1 1 2
+
+let m_1112_1H_2Fa = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+      sum mos_cabs (fun a ->
+        h_one j a Spin.Alfa *. f_two a i j k Spin.Alfa Spin.Alfa ))
+
+    
+let m_1112_1H_2Fb = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+      sum mos_cabs (fun a ->
+        h_one j a Spin.Beta *. f_two a i j k Spin.Alfa Spin.Beta ))
+
+
+let m_1112_2Ha_2Fa =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two j m a m Spin.Alfa Spin.Alfa *.
+      f_two a i j k Spin.Alfa Spin.Alfa 
+    )
+  )
+
+
+let m_1112_2Hb_2Fa =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two j m a m Spin.Alfa Spin.Beta *.
+      f_two a i j k Spin.Alfa Spin.Alfa 
+    )
+  )
+
+let m_1112_2Ha_2Fb =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two j m a m Spin.Alfa Spin.Alfa *.
+      f_two a i j k Spin.Alfa Spin.Beta 
+    )
+  )
+
+
+let m_1112_2Hb_2Fb =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two j m a m Spin.Alfa Spin.Beta *.
+      f_two a i j k Spin.Alfa Spin.Beta 
+    )
+  )
+
+
+let integral_value 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_j, mos_j' =
+      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)
+      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)
+      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 result = 
+     sum mos_j (fun j -> 
+       m_1112_1H_2Fa.{i,j,k} +. 
+       sum mos_i  (fun m -> m_1112_2Ha_2Fa.{i,j,k,m}) +.
+       sum mos_i' (fun m -> m_1112_2Hb_2Fa.{i,j,k,m})) +.
+     sum mos_j' (fun j ->
+       m_1112_1H_2Fb.{i,j,k} +. 
+       sum mos_i  (fun m -> m_1112_2Hb_2Fb.{i,j,k,m}) +.
+       sum mos_i' (fun m -> m_1112_2Ha_2Fb.{i,j,k,m}))
+    in
+    
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+   
+   
+
+let _ = check 0 integral_value 1 1 1 2
+
+let _ = png_image "1_1_21.png"
+
+let integral_value 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_j, mos_j' =
+      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)
+      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)
+      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 kj, mos_b kj
+      | Spin.Beta -> mos_b kj, mos_a kj
+    in
+
+    let result = 
+       let s' = Spin.other s in
+       sum mos_cabs (fun a -> 
+       sum mos_j    (fun j -> 
+               (h_two i j k a s s ) *.
+               (f_one a j s  +. 
+                  sum mos_i  (fun m -> f_two m a m j s  s) +.
+                  sum mos_i' (fun m -> f_two m a m j s' s) ) 
+             ) 
+       +.
+       sum mos_j' (fun j -> 
+               (h_two i j k a s s') *.
+               (f_one a j s' +. 
+                  sum mos_i  (fun m -> f_two m a m j s  s') +.
+                  sum mos_i' (fun m -> f_two m a m j s' s') ) 
+       )) 
+    in
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 1 1 2 1
+
+let m_1121_2Ha_1F = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+      sum mos_cabs (fun a ->
+        h_two i j k a Spin.Alfa Spin.Alfa *. f_one a j Spin.Alfa))
+
+    
+let m_1121_2Hb_1F = 
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+      sum mos_cabs (fun a ->
+        h_two i j k a Spin.Alfa Spin.Beta *. f_one a j Spin.Beta))
+
+
+let m_1121_2Ha_2Fa =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two i j k a Spin.Alfa Spin.Alfa *.
+      f_two m a m j Spin.Alfa Spin.Alfa 
+    )
+  )
+
+
+let m_1121_2Hb_2Fa =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two i j k a Spin.Alfa Spin.Beta *.
+      f_two m a m j Spin.Alfa Spin.Alfa 
+    )
+  )
+
+let m_1121_2Ha_2Fb =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two i j k a Spin.Alfa Spin.Alfa *.
+      f_two m a m j Spin.Alfa Spin.Beta 
+    )
+  )
+
+
+let m_1121_2Hb_2Fb =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two i j k a Spin.Alfa Spin.Beta *.
+      f_two m a m j Spin.Alfa Spin.Beta 
+    )
+  )
+
+
+let integral_value 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_j, mos_j' =
+      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)
+      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)
+      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 kj, mos_b kj
+      | Spin.Beta -> mos_b kj, mos_a kj
+    in
+   
+    let result = 
+     sum mos_j (fun j -> 
+       m_1121_2Ha_1F.{i,j,k} +. 
+       sum mos_i  (fun m -> m_1121_2Ha_2Fa.{i,j,k,m}) +.
+       sum mos_i' (fun m -> m_1121_2Ha_2Fb.{i,j,k,m})) +.
+     sum mos_j' (fun j ->
+       m_1121_2Hb_1F.{i,j,k} +. 
+       sum mos_i  (fun m -> m_1121_2Hb_2Fb.{i,j,k,m}) +.
+       sum mos_i' (fun m -> m_1121_2Hb_2Fa.{i,j,k,m}))
+    in
+    
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+   
+   
+
+let _ = check 0 integral_value 1 1 2 1
+
+let _ = png_image "1_1_22.png"
+
+(*
+let ki = det_I.(6)
+let kj = det_I.(58)
+
+let _ = 
+    Format.printf "|I> -> |J> : %a |\n@." Excitation.pp (Excitation.of_det ki kj) ;
+generate_alphas ki kj 1 1 2 2
+|> Array.of_list
+|> Array.mapi (fun kk alpha -> 
+
+(*
+let _ = Determinant.to_lists ki
+let _ = Determinant.to_lists alpha
+let _ = Determinant.to_lists kj
+
+   
+let _ = 
+    Format.printf "|I> -> |J> : %a\n@." Excitation.pp (Excitation.of_det ki kj);
+    Format.printf "|I> -> |a> : %a\n@." Excitation.pp (Excitation.of_det ki alpha);
+    Format.printf "|a> -> |J> : %a\n@." Excitation.pp (Excitation.of_det alpha kj)
+*)
+
+    let integral_value ki kj = 
+    
+      let exc0 = Array.init (aux_num+1) (fun _ -> [|"-";"-"|]) in
+      let exc1 = Array.init (aux_num+1) (fun _ -> [|"-";"-"|]) in 
+      let exc2 = Array.init (aux_num+1) (fun _ -> [|"-";"-"|]) in 
+ 
+      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 spin = function
+        | Spin.Alfa -> 0
+        | _ -> 1
+      in
+        exc0.(i).(spin s ) <- "i" 
+      ; exc0.(k).(spin s ) <- "k" 
+      ;
+      let s0 = s in
+    
+      let i, j, k, l, s, s', p1 =
+          match Excitation.of_det ki alpha with
+          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
+            hole, hole', particle, particle', spin, spin', phase
+          | _ -> assert false
+      in
+      
+      if exc0.(i).(spin s ) =  "-" then exc0.(i).(spin s ) <- "n";
+      if exc0.(j).(spin s') =  "-" then exc0.(j).(spin s') <- "n";
+      if exc0.(k).(spin s ) =  "-" then exc0.(k).(spin s ) <- if k > mo_num then "a" else "m";
+      if exc0.(l).(spin s') =  "-" then exc0.(l).(spin s') <- if l > mo_num then "a" else "m";
+      
+      let string_h = 
+          Printf.sprintf "h_two %s %s %s %s %s %s *.  " 
+            exc0.(i).(spin s )
+            exc0.(j).(spin s')
+            exc0.(k).(spin s )
+            exc0.(l).(spin s')
+            (if s  = s0 then "s " else "s'")
+            (if s' = s0 then "s " else "s'")
+            (*
+            (if exc0.(i).(spin s ) = "i" || exc0.(k).(spin s ) = "k" then "s  " else 
+             if exc0.(i).(spin s ) = "j" || exc0.(k).(spin s ) = "l" then "s' " else "s''")
+            (if exc0.(j).(spin s') = "i" || exc0.(l).(spin s') = "k" then "s  " else 
+             if exc0.(j).(spin s') = "j" || exc0.(l).(spin s') = "l" then "s' " else "s''")
+             *)
+      in 
+      
+      let i, j, k, l, s, s', p2 =
+          match Excitation.of_det alpha 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 string_f = 
+          Printf.sprintf "f_two %s %s %s %s %s %s" 
+            exc0.(i).(spin s )
+            exc0.(j).(spin s')
+            exc0.(k).(spin s )
+            exc0.(l).(spin s')
+            (*
+            (if exc0.(i).(spin s ) = "i" || exc0.(k).(spin s ) = "k" then "s  " else 
+             if exc0.(i).(spin s ) = "j" || exc0.(k).(spin s ) = "l" then "s' " else 
+             if exc0.(i).(spin s ) = "n" || exc0.(k).(spin s ) = "n" then "s''" else
+             if s  = s0 then "s" else "s'")
+            (if exc0.(j).(spin s') = "i" || exc0.(l).(spin s') = "k" then "s  " else 
+             if exc0.(j).(spin s') = "j" || exc0.(l).(spin s') = "l" then "s' " else 
+             if exc0.(j).(spin s') = "n" || exc0.(l).(spin s') = "n" then "s''" else
+             if s' = s0 then "s" else "s'")
+             *)
+            (if s  = s0 then "s " else "s'")
+            (if s' = s0 then "s " else "s'")
+      in
+      
+    Format.printf "|I> -> |a> : %a | %s\n@." Excitation.pp (Excitation.of_det ki alpha) string_h ;
+    Format.printf "|a> -> |J> : %a | %s\n@." Excitation.pp (Excitation.of_det alpha kj) string_f ;
+    
+    let sign = 
+        if Phase.add p1 p2 = phase then "+. " else "-. "
+    in
+sign ^ string_h ^ string_f
+(*
+        
+        let mos, mos' =
+          match s with
+          | Spin.Alfa -> mos_a ki, mos_b ki
+          | Spin.Beta -> mos_b ki, mos_a ki
+        in
+
+        let result = 
+           let k=4 and l=5 in
+           let i=1 and j=2 in
+           let n=2 and a=7 in
+           let s   = Spin.Alfa
+           and s'  = Spin.Alfa
+           and s'' = Spin.Beta
+           in
+           (* 
+           h_two j n l a s' s'' *. f_two a i n k s s''
+           *)
+           h_two i n k a s s'' *. f_two j a l n s s''
+        in
+
+        match phase with
+        | Phase.Pos -> result
+        | Phase.Neg -> -. result
+        *)
+
+in
+integral_value ki kj
+(*
+let a = (compute_HaaF ki [alpha] kj)
+and b = (integral_value ki kj)
+in
+if kk = 31 then
+    Format.printf "%6d %e %e@.@." (kk) a b
+    *)
+)
+|> Array.to_list
+|> List.sort_uniq compare
+
+let _ = compute_HaaF ki (generate_alphas ki kj 1 1 2 2) kj
+
+let _ = integral_value ki kj
+
+*)
+
+let integral_value 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
+    
+   (* MOs unoccupied in both I and J *)
+   let mos_virt_a, mos_virt_b = 
+     Array.init mo_num (fun i -> Some (i+1)) , 
+     Array.init mo_num (fun i -> Some (i+1)) 
+   in
+   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
+   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
+   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
+   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
+   
+   let mos_virt_a, mos_virt_b = 
+     Array.to_list mos_virt_a |> Util.list_some,
+     Array.to_list mos_virt_b |> Util.list_some
+   in
+
+   let mos_virt, mos_virt' =
+      match s with
+      | Spin.Alfa -> mos_virt_a, mos_virt_b
+      | Spin.Beta -> mos_virt_b, mos_virt_a
+   in
+    
+   let mos_j, mos_j' =
+      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)
+      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)
+      in
+      match s with
+      | Spin.Alfa -> alfa, beta
+      | Spin.Beta -> beta, alfa
+   in
+   let result = 
+       let s' = Spin.other s in
+        sum mos_cabs (fun a -> 
+        sum mos_j    (fun n -> 
+        sum mos_virt (fun m -> 
+             h_two n i a m s  s  *.  f_two m a k n s  s 
+        ))) +.
+        sum mos_cabs (fun a' -> 
+        sum mos_j'   (fun n' -> 
+        sum mos_virt (fun m -> 
+             h_two i n' m a' s  s' *.  f_two m a' k n' s  s'
+        ))) +.
+        sum mos_cabs  (fun a -> 
+        sum mos_j'    (fun n' -> 
+        sum mos_virt' (fun m' -> 
+             h_two i n' a m' s  s' *.  f_two a m' k n' s  s'
+        ))) -. 
+        sum mos_cabs (fun a -> 
+        sum mos_j    (fun n -> 
+        sum mos_j    (fun l -> if l <=n then 0. else
+             h_two n l a k s  s  *.  f_two i a l n s  s 
+        ))) -.  
+        sum mos_cabs (fun a' -> 
+        sum mos_j'   (fun n' -> 
+        sum mos_j    (fun l -> 
+             h_two l n' k a' s s'  *.  f_two i a' l n' s s'
+        ))) 
+    in
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+let _ = check 0 integral_value 1 1 2 2
+
+
+
+let m_1122_va =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two j i a m Spin.Alfa Spin.Alfa *.
+      f_two m a k j Spin.Alfa Spin.Alfa 
+    )
+  )
+
+let m_1122_v2 =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two i j m a Spin.Alfa Spin.Beta *.
+      f_two m a k j Spin.Alfa Spin.Beta 
+    )
+  )
+
+let m_1122_v3 =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k m ->
+    sum mos_cabs (fun a ->
+      h_two i j a m Spin.Alfa Spin.Beta *.
+      f_two a m k j Spin.Alfa Spin.Beta 
+    )
+  )
+
+let m_1122_oa =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+    sum mos_cabs (fun a ->
+      h_two j l a k Spin.Alfa Spin.Alfa *.
+      f_two i a l j Spin.Alfa Spin.Alfa 
+    )
+  )
+
+let m_1122_o =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+    sum mos_cabs (fun a ->
+      h_two l j k a Spin.Alfa Spin.Beta *.
+      f_two i a l j Spin.Alfa Spin.Beta 
+    )
+  )
+
+
+let integral_value 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
+    
+   (* MOs unoccupied in both I and J *)
+   let mos_virt_a, mos_virt_b = 
+     Array.init mo_num (fun i -> Some (i+1)) , 
+     Array.init mo_num (fun i -> Some (i+1)) 
+   in
+   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
+   List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
+   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
+   List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
+   
+   let mos_virt_a, mos_virt_b = 
+     Array.to_list mos_virt_a |> Util.list_some,
+     Array.to_list mos_virt_b |> Util.list_some
+   in
+
+   let mos_virt, mos_virt' =
+      match s with
+      | Spin.Alfa -> mos_virt_a, mos_virt_b
+      | Spin.Beta -> mos_virt_b, mos_virt_a
+   in
+    
+   let mos_j, mos_j' =
+      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)
+      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)
+      in
+      match s with
+      | Spin.Alfa -> alfa, beta
+      | Spin.Beta -> beta, alfa
+   in
+   
+    let result = 
+        sum mos_j    (fun j -> 
+        sum mos_virt (fun m -> 
+          m_1122_va.{i,j,k,m}
+        )) +.
+        sum mos_j'   (fun j -> 
+        sum mos_virt (fun m -> 
+          m_1122_v2.{i,j,k,m}
+        ))   +.
+        sum mos_j'    (fun j -> 
+        sum mos_virt' (fun m -> 
+          m_1122_v3.{i,j,k,m}
+        ))   -.
+        sum mos_j    (fun j -> 
+        sum mos_j    (fun l -> if l <= j then 0. else
+          m_1122_oa.{i,j,k,l}
+        )) -. 
+        sum mos_j'   (fun j -> 
+        sum mos_j    (fun l -> 
+          m_1122_o.{i,j,k,l}
+        )) 
+    in
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 1 1 2 2
+
+
+
+
+let _ = png_image "1_2_22.png"
+
+let integral_value ki kj = 
+    let h, p, s, phase =
+          match Excitation.of_det ki kj with
+          | Excitation.(Single (phase, { hole ; particle ; spin })) ->
+            hole, particle, spin, phase
+          | _ -> assert false
+    in
+
+    let mos, mos' =
+      match s with
+      | Spin.Alfa -> mos_a ki, mos_b ki
+      | Spin.Beta -> mos_b ki, mos_a ki
+    in
+   
+   
+
+    let result = 
+       (* Alpha-Beta *)
+       let s' = Spin.other s in
+       sum mos_cabs (fun b -> 
+       sum mos_cabs (fun a -> 
+       sum mos'     (fun j -> h_two h j a b s s'  *. f_two a b p j s s')
+       ))
+       +.
+       (* Alpha-Alpha / Beta-Beta *)
+       sum mos_cabs (fun b -> 
+       sum mos_cabs (fun a -> if b >= a then 0. else
+       sum mos  (fun j -> h_two h j a b s s  *. f_two a b p j s s)
+       )) 
+    in
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 1 2 2 2
+
+let m_1222a =
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+    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 k j Spin.Alfa Spin.Alfa 
+    )
+    )
+  )
+
+let m_1222 =
+  array_3_init mo_num mo_num mo_num (fun i j k ->
+    sum mos_cabs (fun a ->
+    sum mos_cabs (fun b ->
+      h_two i j a b Spin.Alfa Spin.Beta *.
+      f_two a b k j Spin.Alfa Spin.Beta 
+    )
+    )
+  )
+  
+
+let integral_value 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
+    
+   (* MOs unoccupied in both I and J *)
+    let mos_j, mos_j' =
+      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)
+      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)
+      in
+      match s with
+      | Spin.Alfa -> alfa, beta
+      | Spin.Beta -> beta, alfa
+    in
+
+    let result = 
+        sum mos_j    (fun j ->  m_1222a.{i,j,k} ) +.
+        sum mos_j'   (fun j ->  m_1222.{i,j,k} )
+    in
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 1 2 2 2
+
+
+
+
+let _ = png_image "2_1_12.png"
+
+let integral_value 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, mos' =
+          match s with
+          | Spin.Alfa -> mos_a ki, mos_b ki
+          | Spin.Beta -> mos_b ki, mos_a ki
+        in
+        
+        let mos2, mos2' =
+          match s' with
+          | Spin.Alfa -> mos_a ki, mos_b ki
+          | Spin.Beta -> mos_b ki, mos_a ki
+        in
+
+        let result = 
+       sum mos_cabs (fun a -> 
+           let s'' = Spin.other s in
+             ( h_one i a s +.
+               sum mos  (fun n -> h_two i n a n s s  ) +.
+               sum mos' (fun n -> h_two i n a n s s'') 
+             ) *. f_two a j k l s s'
+             ) +.
+       sum mos_cabs (fun a -> 
+           let s'' = Spin.other s' in
+             ( h_one j a s' +.
+               sum mos2 (fun n -> h_two j n a n s' s' ) +.
+               sum mos2'  (fun n -> h_two j n a n s' s'') 
+             ) *. f_two i a k l s s'
+             ) 
+        in
+
+        match phase with
+        | Phase.Pos -> result
+        | Phase.Neg -> -. result
+
+
+let _ = check 0 integral_value 2 1 1 2
+
+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 i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Alfa +.
+        h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Alfa )
+    )
+
+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 i a Spin.Alfa *. f_two a j k l Spin.Alfa Spin.Beta +.
+        h_one j a Spin.Alfa *. f_two i a k l Spin.Alfa Spin.Beta)
+    )
+
+let m_2112_2Ha_2Fa =
+  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
+    sum mos_cabs (fun a ->
+        h_two i n a n Spin.Alfa Spin.Alfa *.
+        f_two a j k l Spin.Alfa Spin.Alfa )
+    )
+    
+let m_2112_2Hb_2Fa =
+  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
+    sum mos_cabs (fun a ->
+        h_two i n a n Spin.Alfa Spin.Beta *.
+        f_two a j k l Spin.Alfa Spin.Alfa )
+    )
+    
+let m_2112_2Ha_2Fb =
+  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
+    sum mos_cabs (fun a ->
+        h_two i n a n Spin.Alfa Spin.Alfa *.
+        f_two a j k l Spin.Alfa Spin.Beta )
+    )
+    
+let m_2112_2Hb_2Fb =
+  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
+    sum mos_cabs (fun a ->
+        h_two i n a n Spin.Alfa Spin.Beta *.
+        f_two a j k l Spin.Alfa Spin.Beta )
+    )
+    
+let integral_value 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, mos' =
+       match s with
+       | Spin.Alfa -> mos_a ki, mos_b ki
+       | Spin.Beta -> mos_b ki, mos_a ki
+     in
+     
+     let result = 
+       if s = s' then
+         m_2112_1H_2Fa.{i,j,k,l} +. 
+         sum mos  (fun n ->
+             m_2112_2Ha_2Fa.{i,j,k,l,n} +. m_2112_2Ha_2Fa.{j,i,l,k,n}
+         ) +.
+         sum mos' (fun n ->
+             m_2112_2Hb_2Fa.{i,j,k,l,n} +. m_2112_2Hb_2Fa.{j,i,l,k,n}
+         ) 
+       else
+         m_2112_1H_2Fb.{i,j,k,l} +.
+         sum mos  (fun n ->
+             m_2112_2Ha_2Fb.{i,j,k,l,n} +. m_2112_2Hb_2Fb.{j,i,l,k,n}
+         ) +.
+         sum mos' (fun n ->
+             m_2112_2Hb_2Fb.{i,j,k,l,n} +. m_2112_2Ha_2Fb.{j,i,l,k,n}
+         ) 
+     in
+
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 2 1 1 2
+
+
+
+
+let _ = png_image "2_1_21.png"
+
+let integral_value 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 result = 
+        
+        let mos, mos', s'' =
+          match s with
+          | Spin.Alfa -> mos_a kj, mos_b kj, Spin.Beta
+          | Spin.Beta -> mos_b kj, mos_a kj, Spin.Alfa
+        in
+        sum mos_cabs (fun a -> 
+             h_two i j a l s s' *.
+             ( f_one a k s +.
+               sum mos  (fun n -> f_two a n k n s s) +.
+               sum mos' (fun n -> f_two a n k n s s'') 
+             ) ) +.
+
+        let mos, mos', s'' =
+          match s' with
+          | Spin.Alfa -> mos_a kj, mos_b kj, Spin.Beta
+          | Spin.Beta -> mos_b kj, mos_a kj, Spin.Alfa
+        in
+        sum mos_cabs (fun a -> 
+             h_two j i a k s' s *.
+             ( f_one a l s' +.
+               sum mos  (fun n -> f_two a n l n s' s') +.
+               sum mos' (fun n -> f_two a n l n s' s'') 
+             ) 
+             ) 
+        in
+
+        match phase with
+        | Phase.Pos -> result
+        | Phase.Neg -> -. result
+
+
+let _ = check 0 integral_value 2 1 2 1
+
+let m_2121_2Ha_1F =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+    sum mos_cabs (fun a ->
+        h_two i j a l Spin.Alfa Spin.Alfa *. f_one a k Spin.Alfa +.
+        h_two i j k a Spin.Alfa Spin.Alfa *. f_one a l Spin.Alfa)
+    )
+
+let m_2121_2Hb_1F =
+  array_4_init mo_num mo_num mo_num mo_num (fun i j k l ->
+    sum mos_cabs (fun a ->
+        h_two i j a l Spin.Alfa Spin.Beta *. f_one a k Spin.Alfa +.
+        h_two i j k a Spin.Alfa Spin.Beta *. f_one a l Spin.Alfa)
+    )
+
+let m_2121_2Ha_2Fa =
+  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
+    sum mos_cabs (fun a ->
+        h_two i j a l Spin.Alfa Spin.Alfa *.
+        f_two a n k n Spin.Alfa Spin.Alfa )
+    )
+    
+let m_2121_2Hb_2Fa =
+  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
+    sum mos_cabs (fun a ->
+        h_two i j a l Spin.Alfa Spin.Beta *.
+        f_two a n k n Spin.Alfa Spin.Alfa )
+    )
+    
+let m_2121_2Ha_2Fb =
+  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
+    sum mos_cabs (fun a ->
+        h_two i j a l Spin.Alfa Spin.Alfa *.
+        f_two a n k n Spin.Alfa Spin.Beta )
+    )
+    
+let m_2121_2Hb_2Fb =
+  array_5_init mo_num mo_num mo_num mo_num mo_num (fun i j k l n ->
+    sum mos_cabs (fun a ->
+        h_two i j a l Spin.Alfa Spin.Beta *.
+        f_two a n k n Spin.Alfa Spin.Beta )
+    )
+    
+let integral_value 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, mos' =
+       match s with
+       | Spin.Alfa -> mos_a kj, mos_b kj
+       | Spin.Beta -> mos_b kj, mos_a kj
+     in
+     
+     let result = 
+       if s = s' then
+         m_2121_2Ha_1F.{i,j,k,l} +. 
+         sum mos  (fun n ->
+             m_2121_2Ha_2Fa.{i,j,k,l,n} +. m_2121_2Ha_2Fa.{j,i,l,k,n}
+         ) +.
+         sum mos' (fun n ->
+             m_2121_2Ha_2Fb.{i,j,k,l,n} +. m_2121_2Ha_2Fb.{j,i,l,k,n}
+         ) 
+       else
+         m_2121_2Hb_1F.{i,j,k,l} +.
+         sum mos  (fun n ->
+             m_2121_2Hb_2Fa.{i,j,k,l,n} +. m_2121_2Hb_2Fb.{j,i,l,k,n}
+         ) +.
+         sum mos' (fun n ->
+             m_2121_2Hb_2Fb.{i,j,k,l,n} +. m_2121_2Hb_2Fa.{j,i,l,k,n}
+         ) 
+     in
+
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 2 1 2 1
+
+
+
+
+let _ = png_image "2_1_22.png"
+
+let integral_value 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_virt_a, mos_virt_b = 
+       Array.init mo_num (fun i -> Some (i+1)) , 
+       Array.init mo_num (fun i -> Some (i+1)) 
+    in
+    List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a ki);
+    List.iter (fun i -> mos_virt_a.(i-1) <- None) (mos_a kj);
+    List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b ki);
+    List.iter (fun i -> mos_virt_b.(i-1) <- None) (mos_b kj);
+
+    let mos_virt_a, mos_virt_b = 
+      Array.to_list mos_virt_a |> Util.list_some,
+      Array.to_list mos_virt_b |> Util.list_some
+    in
+   
+    let result = 
+        
+        let mos_virt, mos_virt' =
+          match s with
+          | Spin.Alfa -> mos_virt_a, mos_virt_b
+          | Spin.Beta -> mos_virt_b, mos_virt_a
+        in
+        
+        let mos, mos' =
+          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)
+          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)
+          in
+          match s with
+          | Spin.Alfa -> alfa, beta
+          | Spin.Beta -> beta, alfa
+        in
+
+        if s = s' then            
+            let s'' = Spin.other s' in
+            sum mos_cabs (fun a -> 
+            sum mos  (fun n ->
+                 h_two i n a k s  s  *.  f_two j a n l s  s 
+              +. h_two i n a l s  s  *.  f_two j a k n s  s 
+              -. h_two j n a k s  s  *.  f_two i a n l s  s   
+              -. h_two j n a l s  s  *.  f_two i a k n s  s
+              )
+            +. sum mos_virt (fun m ->
+                -. h_two i j a m s  s  *.  f_two m a k l s  s )
+            +. sum mos' (fun n ->
+                 h_two i n k a s  s'' *.  f_two j a l n s  s''
+              +. h_two j n l a s  s'' *.  f_two i a k n s  s''
+              -. h_two i n l a s  s'' *.  f_two j a k n s  s''
+              -. h_two j n k a s  s'' *.  f_two i a l n s  s''
+              )
+            )
+        else
+            sum mos_cabs (fun a ->
+            sum mos_virt' (fun m ->
+                  h_two i j a m s  s' *.  f_two a m k l s  s' ) +.
+            sum mos_virt  (fun m ->
+                  h_two i j m a s  s' *.  f_two m a k l s  s' ) +.
+            sum mos  (fun n ->
+                  h_two n i a k s  s  *.  f_two a j n l s  s'
+               +. h_two n j a l s  s' *.  f_two i a k n s  s 
+               -. h_two n j k a s  s' *.  f_two i a n l s  s'
+               ) +.
+            sum mos' (fun n -> if n >= j then 0. else 
+                  h_two i n k a s  s' *.  f_two j a l n s' s'
+               +. h_two n j a l s' s' *.  f_two i a k n s  s' ) +. 
+            sum mos' (fun n -> if n <= j then 0. else 
+               -. h_two i n k a s  s' *.  f_two j a n l s' s' 
+               -. h_two j n a l s' s' *.  f_two i a k n s  s' ) +.
+            sum mos' (fun n -> 
+               -. h_two i n a l s  s' *.  f_two a j k n s  s'
+               )
+            )
+    in
+
+    match phase with
+        | Phase.Pos -> result
+        | Phase.Neg -> -. result
+
+
+let _ = check 0 integral_value 2 1 2 2
+
+let _ = png_image "2_2_22.png"
+
+let integral_value ki kj = 
+    let h, h', p, p', s, s', phase =
+          match Excitation.of_det ki kj with
+          | Excitation.(Double (phase,
+              { hole=h ; particle=p ; spin=s },
+              { hole=h'; particle=p'; spin=s'}) )  -> h, h', p, p', s, s', phase
+          | _ -> assert false
+    in
+
+    let result = 
+       if s <> s' then (* Alpha-Beta *)
+           sum mos_cabs (fun b -> 
+           sum mos_cabs (fun a -> 
+               h_two h h' a b s s'  *. f_two a b p p' s s'
+           )) 
+       else (* Alpha-Alpha / Beta-Beta *)
+           sum mos_cabs (fun b -> 
+           sum mos_cabs (fun a ->  if b >= a then 0. else
+               h_two h h' a b s s'  *. f_two a b p p' s s'
+           )) 
+    in
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 2 2 2 2
+
+let _ = png_image "3_1_22.png"
+
+let integral_value 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 i j a k s1 s2 *.  f_two m a l n s3 s3 
+             +. h_two i j a n s1 s2 *.  f_two m a k l s3 s2 
+             +. h_two i m a l s1 s3 *.  f_two j a k n s2 s3 
+             +. h_two j m a k s2 s3 *.  f_two i a l n s1 s3 
+             +. h_two j m a n s2 s3 *.  f_two i a k l s1 s2 
+             -. h_two i j a l s1 s2 *.  f_two m a k n s3 s3 
+             -. h_two i m a k s1 s3 *.  f_two j a l n s2 s3 
+             -. h_two i m a n s1 s3 *.  f_two j a k l s2 s2 
+             -. h_two j m a l s2 s3 *.  f_two i a k n s1 s3 )
+       | Alfa, Alfa, Beta
+       | Beta, Beta, Alfa -> 
+           sum mos_cabs (fun a -> 
+                h_two i j a l s1 s2 *.  f_two a m k n s1 s3
+             +. h_two i m k a s1 s3 *.  f_two j a l n s2 s3
+             +. h_two j m a n s2 s3 *.  f_two i a k l s1 s2 
+             +. h_two j m l a s2 s3 *.  f_two i a k n s1 s3
+             -. h_two i j a k s1 s2 *.  f_two a m l n s1 s3
+             -. h_two i m a n s1 s3 *.  f_two j a k l s2 s2 
+             -. h_two i m l a s1 s3 *.  f_two j a k n s2 s3
+             -. h_two j m k a s2 s3 *.  f_two i a l n s1 s3
+           ) 
+       | Alfa, Beta, Beta
+       | Beta, Alfa, Alfa -> 
+           sum mos_cabs (fun a -> 
+                h_two i j a l s1 s2 *.  f_two a m k n s1 s3 
+             +. h_two i m a n s1 s3 *.  f_two a j k l s1 s2 
+             +. h_two i m k a s1 s3 *.  f_two j a l n s2 s3 
+             +. h_two j m a n s2 s3 *.  f_two i a k l s1 s2 
+             -. h_two i j a n s1 s2 *.  f_two a m k l s1 s2 
+             -. h_two i j k a s1 s2 *.  f_two m a l n s2 s3 
+             -. h_two i m a l s1 s3 *.  f_two a j k n s1 s3 
+             -. h_two j m a l s2 s3 *.  f_two i a k n s1 s3 
+           ) 
+       | Beta, Alfa, Beta
+       | Alfa, Beta, Alfa ->  0. (*TODO *)
+    in
+    match phase with
+    | Phase.Pos -> result
+    | Phase.Neg -> -. result
+     
+
+let _ = check 0 integral_value 3 1 2 2
+
+let ki = det_I.(129)
+let kj = det_I.(349)
+
+let alpha_to_string alpha = 
+   
+      let exc0 = Array.init (aux_num+1) (fun _ -> [|"-";"-"|]) in
+ 
+      let i, j, m, k, l, n,  s, s', s'',  phase =
+          match Excitation.of_det ki kj with
+          | Excitation.(Triple (phase,
+              { hole ; particle ; spin },
+              {hole=hole' ; particle=particle' ; spin=spin' }, 
+              {hole=hole''; particle=particle''; spin=spin''}  )) ->
+            hole, hole', hole'', particle, particle', particle'', spin, spin', spin'', phase
+          | _ -> assert false
+      in
+      let spin = function
+        | Spin.Alfa -> 0
+        | _ -> 1
+      in
+        exc0.(i).(spin s  ) <- "i" 
+      ; exc0.(j).(spin s' ) <- "j" 
+      ; exc0.(k).(spin s  ) <- "k" 
+      ; exc0.(l).(spin s' ) <- "l" 
+      ; exc0.(m).(spin s'') <- "m" 
+      ; exc0.(n).(spin s'') <- "n" 
+      ;
+      let s0, s0', s0'' = s, s', s'' in
+    
+      let i, j, k, l, s, s', p1 =
+          match Excitation.of_det ki alpha with
+          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
+            hole, hole', particle, particle', spin, spin', phase
+          | _ -> assert false
+      in
+      
+      if exc0.(i).(spin s ) =  "-" then exc0.(i).(spin s ) <- "p";
+      if exc0.(j).(spin s') =  "-" then exc0.(j).(spin s') <- "p";
+      if exc0.(k).(spin s ) =  "-" then exc0.(k).(spin s ) <- if k > mo_num then "a" else "q";
+      if exc0.(l).(spin s') =  "-" then exc0.(l).(spin s') <- if l > mo_num then "a" else "q";
+      
+      let string_h = 
+          Printf.sprintf "h_two %s %s %s %s %s %s *.  " 
+            exc0.(i).(spin s )
+            exc0.(j).(spin s')
+            exc0.(k).(spin s )
+            exc0.(l).(spin s')
+            (if s  = s0  then "s " else if s = s0'  then "s'"  else "s''")
+            (if s' = s0' then "s'" else if s = s0'' then "s''" else "s"  )
+      in 
+      
+      let i, j, k, l, s, s', p2 =
+          match Excitation.of_det alpha 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 string_f = 
+          Printf.sprintf "f_two %s %s %s %s %s %s" 
+            exc0.(i).(spin s )
+            exc0.(j).(spin s')
+            exc0.(k).(spin s )
+            exc0.(l).(spin s')
+            (if s  = s0  then "s " else if s = s0'  then "s'"  else "s''")
+            (if s' = s0' then "s'" else if s = s0'' then "s''" else "s"  )
+      in
+        (*
+        Format.printf "|I> -> |a> : %a | %s\n@." Excitation.pp (Excitation.of_det ki alpha) string_h ;
+        Format.printf "|a> -> |J> : %a | %s\n@." Excitation.pp (Excitation.of_det alpha kj) string_f ;
+        *)
+
+    
+        let sign = 
+            if Phase.add p1 p2 = phase then "+." else "-."
+        in
+        sign ^ string_h ^ string_f
+    
+
+let alpha_debug alpha = 
+   
+      let i, j, k, l, s, s', p1 =
+          match Excitation.of_det ki alpha with
+          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
+            hole, hole', particle, particle', spin, spin', phase
+          | _ -> assert false
+      in
+      Printf.printf "%d %d %d %d " i j k l;
+      
+      let i, j, k, l, s, s', p2 =
+          match Excitation.of_det alpha kj with
+          | Excitation.(Double (phase, { hole ; particle ; spin }, {hole=hole' ; particle=particle' ; spin=spin' })) ->
+            hole, hole', particle, particle', spin, spin', phase
+          | _ -> assert false
+      in
+      (*
+        Format.printf "|I> -> |a> : %a | \n@." Excitation.pp (Excitation.of_det ki alpha) ;
+        Format.printf "|a> -> |J> : %a | \n@." Excitation.pp (Excitation.of_det alpha kj) ;
+        *)
+      Printf.printf "%d %d %d %d \n%!" i j k l
+
+      
+let strings = 
+    Format.printf "|I> -> |J> : %a |\n@." Excitation.pp (Excitation.of_det ki kj) ;
+    generate_alphas ki kj 3 1 2 2
+    |> Array.of_list
+    |> Array.mapi (fun kk alpha -> alpha_to_string alpha)
+    |> Array.to_list
+    |> List.sort_uniq compare
+    |> Array.of_list
+
+let _ = Array.iteri (fun i x -> Printf.printf "%d %s \n%!" i x) strings
+
+let _ =
+  let v =
+    let alphas =
+      generate_alphas ki kj 3 1 2 2
+      (*
+      |> List.filter (fun alpha ->
+        let x =  alpha_to_string alpha in
+           x = strings.(6)
+        )
+        *)
+    in
+    (*
+    List.iter alpha_debug alphas ;
+    Printf.printf "\n%!";
+    *)
+    compute_HaaF ki alphas kj 
+  in
+  let x = (integral_value ki kj) in
+  Printf.printf "%20.8e %20.8e %20.8e\n%!" x v (v-. x)
+
+
+