diff --git a/CI/CI.ml b/CI/CI.ml
index 5ac3280..3afd07c 100644
--- a/CI/CI.ml
+++ b/CI/CI.ml
@@ -22,27 +22,27 @@ let m_H          t = Lazy.force t.m_H
 
 let m_S2         t = Lazy.force t.m_S2
 
-let eigensystem  t = Lazy.force t.eigensystem 
+let eigensystem  t = Lazy.force t.eigensystem
 
 let mo_class     t = DeterminantSpace.mo_class t.det_space
 
-let eigenvectors t = 
+let eigenvectors t =
   let (x,_) = eigensystem t in x
 
-let eigenvalues t = 
+let eigenvalues t =
   let (_,x) = eigensystem t in x
 
-let h_integrals mo_basis = 
+let h_integrals mo_basis =
   let one_e_ints    = MOBasis.one_e_ints mo_basis
   and two_e_ints    = MOBasis.two_e_ints mo_basis
   in
-  ( (fun i j _ -> one_e_ints.{i,j}), 
+  ( (fun i j _ -> one_e_ints.{i,j}),
     (fun i j k l s s' ->
        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) 
+         (ERI.get_phys two_e_ints i j k l) -.
+         (ERI.get_phys two_e_ints i j l k)
     ), None )
 
 
@@ -52,7 +52,7 @@ let h_ij mo_basis ki kj =
     List.map (fun f -> f mo_basis)
       [ h_integrals ]
   in
-  CIMatrixElement.make integrals ki kj 
+  CIMatrixElement.make integrals ki kj
   |> List.hd
 
 
@@ -61,14 +61,14 @@ let h_ij_non_zero mo_basis deg_a deg_b ki kj =
     List.map (fun f -> f mo_basis)
       [ h_integrals ]
   in
-  CIMatrixElement.non_zero integrals deg_a deg_b ki kj 
+  CIMatrixElement.non_zero integrals deg_a deg_b ki kj
   |> List.hd
 
 
 let create_matrix_arbitrary f det_space =
   lazy (
     let ndet = Ds.size det_space in
-    let data = 
+    let data =
       match Ds.determinants det_space with
       | Ds.Arbitrary a -> a
       | _ -> assert false
@@ -81,11 +81,11 @@ let create_matrix_arbitrary f det_space =
 
 
     (** Array of (list of singles, list of doubles) in the beta spin *)
-    let degree_bb = 
+    let degree_bb =
       Array.map (fun det_i ->
           let deg = Spindeterminant.degree det_i in
           let doubles =
-            Array.mapi (fun i det_j -> 
+            Array.mapi (fun i det_j ->
                 let d = deg det_j in
                 if d < 3 then
                   Some (i,d,det_j)
@@ -95,18 +95,18 @@ let create_matrix_arbitrary f det_space =
             |> Array.to_list
             |> Util.list_some
           in
-          let singles = 
+          let singles =
             List.filter (fun (i,d,det_j) -> d < 2) doubles
             |> List.map (fun (i,_,det_j) -> (i,det_j))
           in
-          let doubles = 
+          let doubles =
             List.map (fun (i,_,det_j) -> (i,det_j)) doubles
           in
           (singles, doubles)
         ) det_beta
     in
 
-    let task (i,i_dets) = 
+    let task (i,i_dets) =
       let shift = index_start.(i) in
 
       let result =
@@ -130,11 +130,11 @@ let create_matrix_arbitrary f det_space =
 
           begin
             match degree_a with
-            | 2 -> 
+            | 2 ->
               Array.iteri (fun i' i_b ->
                   try
                     Array.iteri (fun j' j_b ->
-                        if j_b >= i_b then 
+                        if j_b >= i_b then
                           ( if j_b = i_b then
                               ( let i_beta = det_beta.(i_b) in
                                 let ki = Determinant.of_spindeterminants i_alfa i_beta in
@@ -144,12 +144,12 @@ let create_matrix_arbitrary f det_space =
                             raise Exit)
                       ) j_dets
                   with Exit -> ()
-                ) i_dets 
-            | 1 -> 
+                ) i_dets
+            | 1 ->
               Array.iteri (fun i' i_b ->
                   let i_beta = det_beta.(i_b) in
                   let ki = Determinant.of_spindeterminants i_alfa i_beta in
-                  let singles, _ = degree_bb.(i_b) in 
+                  let singles, _ = degree_bb.(i_b) in
                   let rec aux singles j' =
                     match singles with
                     | [] -> ()
@@ -167,13 +167,13 @@ let create_matrix_arbitrary f det_space =
                         |  1 -> if (j' < Array.length j_dets) then (aux [@tailcall]) singles (j'+1)
                         |  _ -> assert false
                       end
-                  in aux singles 0 
-                ) i_dets 
-            | 0 -> 
+                  in aux singles 0
+                ) i_dets
+            | 0 ->
               Array.iteri (fun i' i_b ->
                   let i_beta = det_beta.(i_b) in
                   let ki = Determinant.of_spindeterminants i_alfa i_beta in
-                  let _, doubles = degree_bb.(i_b) in 
+                  let _, doubles = degree_bb.(i_b) in
                   let rec aux doubles j' =
                     match doubles with
                     | [] -> ()
@@ -191,13 +191,13 @@ let create_matrix_arbitrary f det_space =
                         |  1 -> if (j' < Array.length j_dets) then (aux [@tailcall]) doubles (j'+1)
                         |  _ -> assert false
                       end
-                  in aux doubles 0 
-                ) i_dets 
+                  in aux doubles 0
+                ) i_dets
             | _ -> ();
           end
         ) det;
       let r =
-        Array.map (fun l -> 
+        Array.map (fun l ->
             List.rev l
             |> Vector.sparse_of_assoc_list ndet
           ) result
@@ -212,7 +212,7 @@ let create_matrix_arbitrary f det_space =
         Array.init ndet (fun _ -> Vector.sparse_of_assoc_list ndet [])
     in
 
-    let n_det_alfa = 
+    let n_det_alfa =
       Array.length det_alfa
     in
     Array.mapi (fun i i_dets -> i, i_dets) det
@@ -231,10 +231,10 @@ let create_matrix_arbitrary f det_space =
 
 (* Create a matrix using the fact that the determinant space is made of
    the outer product of spindeterminants.  *)
-let create_matrix_spin ?(nmax=2) f det_space = 
+let create_matrix_spin ?(nmax=2) f det_space =
   lazy (
     let ndet = Ds.size det_space in
-    let a, b = 
+    let a, b =
       match Ds.determinants det_space with
       | Ds.Spin (a,b) -> (a,b)
       | _ -> assert false
@@ -243,14 +243,14 @@ let create_matrix_spin ?(nmax=2) f det_space =
 
 
     (** Array of (list of singles, list of doubles, list of triples) in the beta spin *)
-    let degree_bb = 
+    let degree_bb =
       match nmax with
       | 2 -> begin
                 Array.map (fun det_i ->
                     let deg = Spindeterminant.degree det_i in
                     let triples = [] in
                     let doubles =
-                      Array.mapi (fun i det_j -> 
+                      Array.mapi (fun i det_j ->
                           let d = deg det_j in
                           if d < 3 then
                             Some (i,d,det_j)
@@ -260,11 +260,11 @@ let create_matrix_spin ?(nmax=2) f det_space =
                       |> Array.to_list
                       |> Util.list_some
                     in
-                    let singles = 
+                    let singles =
                       List.filter (fun (i,d,det_j) -> d < 2) doubles
                       |> List.map (fun (i,_,det_j) -> (i,det_j))
                     in
-                    let doubles = 
+                    let doubles =
                       List.map (fun (i,_,det_j) -> (i,det_j)) doubles
                     in
                     (singles, doubles, triples)
@@ -273,8 +273,8 @@ let create_matrix_spin ?(nmax=2) f det_space =
       | 3 -> begin
                 Array.map (fun det_i ->
                     let deg = Spindeterminant.degree det_i in
-                    let triples = 
-                      Array.mapi (fun i det_j -> 
+                    let triples =
+                      Array.mapi (fun i det_j ->
                           let d = deg det_j in
                           if d < 4 then
                             Some (i,d,det_j)
@@ -284,20 +284,20 @@ let create_matrix_spin ?(nmax=2) f det_space =
                       |> Array.to_list
                       |> Util.list_some
                     in
-                    let doubles = 
+                    let doubles =
                       List.filter (fun (i,d,det_j) -> d < 3) triples
                     in
-                    let singles = 
+                    let singles =
                       List.filter (fun (i,d,det_j) -> d < 2) doubles
                     in
 
-                    let triples = 
+                    let triples =
                       List.map (fun (i,_,det_j) -> (i,det_j)) triples
                     in
-                    let doubles = 
+                    let doubles =
                       List.map (fun (i,_,det_j) -> (i,det_j)) doubles
                     in
-                    let singles = 
+                    let singles =
                       List.map (fun (i,_,det_j) -> (i,det_j)) singles
                     in
                     (singles, doubles, triples)
@@ -309,7 +309,7 @@ let create_matrix_spin ?(nmax=2) f det_space =
     and b = Array.to_list b
     in
 
-    let task = 
+    let task =
       let result =
         Array.init n_beta (fun _ -> [])
       in
@@ -339,7 +339,7 @@ let create_matrix_spin ?(nmax=2) f det_space =
                               update !i' (ib + !j) 2 0 ki kj;
                               incr i';
                             ) b;
-                        | 1 -> 
+                        | 1 ->
                           let i' = ref 0 in
                           List.iteri (fun ib i_beta ->
                               let ki = Determinant.of_spindeterminants i_alfa i_beta in
@@ -382,7 +382,7 @@ let create_matrix_spin ?(nmax=2) f det_space =
                               update !i' (ib + !j) 3 0 ki kj;
                               incr i';
                             ) b;
-                        | 2 -> 
+                        | 2 ->
                           let i' = ref 0 in
                           List.iteri (fun ib i_beta ->
                               let ki = Determinant.of_spindeterminants i_alfa i_beta in
@@ -422,11 +422,11 @@ let create_matrix_spin ?(nmax=2) f det_space =
                   end
           | _ -> assert false
           end;
-          let r = 
+          let r =
             Array.map (fun l ->
                 List.rev l
                 |> Vector.sparse_of_assoc_list ndet
-              ) result 
+              ) result
           in (i,r)
     in
 
@@ -438,7 +438,7 @@ let create_matrix_spin ?(nmax=2) f det_space =
 
     List.mapi (fun i i_alfa -> i*n_beta, i_alfa) a
     |> Stream.of_list
-    |> Farm.run ~ordered:false ~f:task 
+    |> Farm.run ~ordered:false ~f:task
     |> Stream.iter (fun (k, r) ->
         Array.iteri (fun j r_j -> result.(k+j) <- r_j) r;
         Printf.eprintf "%8d / %8d\r%!" (k+Array.length r) ndet;
@@ -450,10 +450,10 @@ let create_matrix_spin ?(nmax=2) f det_space =
 
 (* Create a matrix using the fact that the determinant space is made of
    the outer product of spindeterminants.  *)
-let create_matrix_spin_computed ?(nmax=2) f det_space = 
+let create_matrix_spin_computed ?(nmax=2) f det_space =
   lazy (
     let ndet = Ds.size det_space in
-    let a, b = 
+    let a, b =
       match Ds.determinants det_space with
       | Ds.Spin (a,b) -> (a,b)
       | _ -> assert false
@@ -579,11 +579,11 @@ let create_matrix_spin_computed ?(nmax=2) f det_space =
           let j1 = ref (2*ndet) in
           let j_a = ref 0 in
           result := fun j ->
-              if  (!j0 < j) && (j <= !j1) then 
+              if  (!j0 < j) && (j <= !j1) then
                   ()
               else
                 begin
-                  if (!j1 < j) && (j <= !j1 + n_beta) then 
+                  if (!j1 < j) && (j <= !j1 + n_beta) then
                     begin
                       incr j_a;
                       j0 := !j1;
@@ -607,7 +607,7 @@ let create_matrix_spin_computed ?(nmax=2) f det_space =
 
 
 
-let make ?(n_states=1) ?(algo=`Direct) det_space = 
+let make ?(n_states=1) ?(algo=`Direct) det_space =
 
   let mo_basis = Ds.mo_basis det_space in
 
@@ -616,23 +616,23 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
   *)
   ignore @@ MOBasis.two_e_ints mo_basis;
 
-  let e_shift = 
-    let d0 = 
+  let e_shift =
+    let d0 =
       Ds.determinant_stream det_space
       |> Stream.next
     in
     h_ij_non_zero mo_basis 0 0 d0 d0
   in
 
-  let m_H = 
+  let m_H =
 
-    let f = 
+    let f =
       match Ds.determinants det_space with
       | Ds.Arbitrary _ -> create_matrix_arbitrary
-      | Ds.Spin      _ -> 
+      | Ds.Spin      _ ->
         if algo = `Direct then
           create_matrix_spin_computed ~nmax:2
-        else 
+        else
           create_matrix_spin  ~nmax:2
     in
     f (fun deg_a deg_b ki kj ->
@@ -640,22 +640,22 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
             h_ij_non_zero mo_basis deg_a deg_b ki kj
          else
             h_ij_non_zero mo_basis 0 0 ki ki -. e_shift
-      ) det_space 
+      ) det_space
   in
 
   let m_S2 =
-    let f = 
+    let f =
       match Ds.determinants det_space with
       | Ds.Arbitrary _ -> create_matrix_arbitrary
       | Ds.Spin      _ -> create_matrix_spin  ~nmax:2
     in
-    f (fun _deg_a _deg_b ki kj -> CIMatrixElement.make_s2 ki kj) det_space 
+    f (fun _deg_a _deg_b ki kj -> CIMatrixElement.make_s2 ki kj) det_space
     (*TODO*)
   in
 
   let eigensystem = lazy (
 
-    let eigensystem_incore () = 
+    let eigensystem_incore () =
       let m_H =
         Lazy.force m_H
       in
@@ -664,15 +664,15 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
           Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i)
         ))
       in
-      let matrix_prod psi = 
-        let result = 
-          Matrix.mm ~transa:`T psi m_H 
+      let matrix_prod psi =
+        let result =
+          Matrix.mm ~transa:`T psi m_H
           |> Matrix.transpose
         in
         Parallel.broadcast (lazy result)
       in
-      let eigenvectors, eigenvalues = 
-        let result = 
+      let eigenvectors, eigenvalues =
+        let result =
           Davidson.make ~threshold:1.e-6  ~n_states diagonal matrix_prod
         in
         Parallel.broadcast (lazy result)
@@ -681,7 +681,7 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
       (Conventions.rephase @@ Util.remove_epsilons eigenvectors), eigenvalues
     in
 
-    let eigensystem_direct () = 
+    let eigensystem_direct () =
       let m_H =
         Lazy.force m_H
       in
@@ -690,15 +690,15 @@ let make ?(n_states=1) ?(algo=`Direct) det_space =
           Vec.init (Matrix.dim1 m_H) (fun i -> Matrix.get m_H i i)
         ))
       in
-      let matrix_prod psi = 
-        let result = 
+      let matrix_prod psi =
+        let result =
           Matrix.parallel_mm ~transa:`T ~transb:`N psi m_H
           |> Matrix.transpose
         in
         Parallel.broadcast (lazy result)
       in
-      let eigenvectors, eigenvalues = 
-        let result = 
+      let eigenvectors, eigenvalues =
+        let result =
           Davidson.make ~threshold:1.e-6  ~n_states diagonal matrix_prod
         in
         Parallel.broadcast (lazy result)
@@ -725,7 +725,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
   let list_particles = Array.of_list list_particles in
 
   let psi0 =
-    let stream = 
+    let stream =
       Ds.determinant_stream det_space
     in
     Array.init (Ds.size det_space) (fun i ->
@@ -768,7 +768,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
       in aux [] (Array.length psi0 - 1)
     in
 
-    let psi_filtered = 
+    let psi_filtered =
       List.map (fun i -> psi0.(i)) psi_filtered_idx
     in
 
@@ -794,7 +794,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     let det_i = fst psi0.(i) in
     let w_alfa = w_alfa det_i in
 
-    let same_spin = 
+    let same_spin =
       List.fold_left (fun accu spin ->
           accu +.
           Array.fold_left (fun accu particle ->
@@ -806,9 +806,9 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                     in
                     if Determinant.is_none alfa then accu else
 
-                      let single = 
+                      let single =
                         if already_generated alfa then 0. else
-                          w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa 
+                          w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa
                       in
 
                       let double =
@@ -821,7 +821,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                                   if hole' = particle' || hole' = particle || hole' <= hole then
                                     accu
                                   else
-                                      let alfa = 
+                                      let alfa =
                                         Determinant.single_excitation
                                           spin  hole' particle' alfa
                                       in
@@ -829,7 +829,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                                          already_generated alfa then
                                         accu
                                       else
-                                        accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa 
+                                        accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa
                                 ) 0. list_holes
                           ) 0. list_particles
                       in
@@ -854,7 +854,7 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                           accu +.
                           Array.fold_left (fun accu hole' ->
                               if hole' = particle' then accu else
-                                let alfa = 
+                                let alfa =
                                   Determinant.double_excitation
                                     Spin.Alfa hole  particle
                                     Spin.Beta hole' particle' det_i
@@ -863,12 +863,12 @@ let second_order_sum { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                                   already_generated alfa then
                                   accu
                                 else
-                                  accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa 
+                                  accu +. w_alfa alfa *. psi_h_alfa_alfa_h_psi alfa
                           ) 0. list_holes
                       ) 0. list_particles
                   in
 
-                  accu +. double_ab 
+                  accu +. double_ab
             ) 0. list_holes
         ) 0. list_particles
     in
@@ -885,7 +885,7 @@ let second_order_sum2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     list_holes  list_particles i_o1_alfa e0 psi0 =
 
   let psi0 =
-    let stream = 
+    let stream =
       Ds.determinant_stream det_space
     in
     Array.init (Ds.size det_space) (fun i ->
@@ -897,7 +897,7 @@ let second_order_sum2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
     Ds.determinants_array det_space
     |> Array.to_list
     |> List.map (fun det_i ->
-        [ Spin.Alfa ; Spin.Beta ] 
+        [ Spin.Alfa ; Spin.Beta ]
         |>  List.map (fun spin ->
             List.map (fun particle ->
               List.map (fun hole ->
@@ -905,7 +905,7 @@ let second_order_sum2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                     List.map (fun particle' ->
                         List.map (fun hole' ->
                           Determinant.double_excitation
-                            spin  hole particle 
+                            spin  hole particle
                             spin  hole' particle' det_i
                         ) list_holes
                     ) list_particles
@@ -914,7 +914,7 @@ let second_order_sum2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
                     List.map (fun particle' ->
                         List.map (fun hole' ->
                           Determinant.double_excitation
-                            spin  hole particle 
+                            spin  hole particle
                             (Spin.other spin)  hole' particle' det_i
                         ) list_holes
                     ) list_particles
@@ -926,7 +926,7 @@ let second_order_sum2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
             |> List.concat
           )
           |> List.concat
-        )  
+        )
         |> List.concat
         |> List.concat
     |> List.filter (fun alfa -> not (Determinant.is_none alfa))
@@ -935,7 +935,7 @@ let second_order_sum2 { det_space ; m_H ; m_S2 ; eigensystem ; n_states }
 
   List.fold_left (fun accu alfa ->
     let alfa_o2 = i_o1_alfa alfa in
-    let a_h_psi = 
+    let a_h_psi =
       Array.fold_left (fun accu (det,ci) -> ci.{1} *. (alfa_o2 det)) 0. psi0
     in
     accu +. (a_h_psi *. a_h_psi) /. (e0 -. (alfa_o2 alfa))
@@ -973,7 +973,35 @@ let is_internal det_space =
         Bitstring.logand neg_active_mask beta = occ_mask
 
 
-let _pt2_en ci =
+let is_external det_space =
+    let mo_class      = DeterminantSpace.mo_class det_space in
+    let numbits = 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 numbits j)
+        ) (Bitstring.zero numbits) l
+    in
+    let inactive_mask = m (MOClass.inactive_mos mo_class) in
+    fun a ->
+      let alfa =
+        Determinant.alfa a
+        |> Spindeterminant.bitstring
+      in
+      let n_a = 
+	Bitstring.(popcount @@ logand inactive_mask alfa) 
+      in
+      match n_a with
+      | 0 | 1 | 2 -> 
+	let beta =
+	  Determinant.beta a
+	  |> Spindeterminant.bitstring
+	in
+	n_a + Bitstring.(popcount @@ logand inactive_mask beta) < 3 
+      | _ ->   false
+
+
+let pt2_en ci =
 
   let mo_basis  = Ds.mo_basis ci.det_space in
   let psi0, e0 = Parallel.broadcast ci.eigensystem in
@@ -982,14 +1010,14 @@ let _pt2_en ci =
 
   let w_alfa det_i =
     let one_e, two_e, _ = h_integrals mo_basis in
-    let fock_diag_alfa, fock_diag_beta = 
+    let fock_diag_alfa, fock_diag_beta =
       Ds.fock_diag ci.det_space det_i
     in
     let h_aa alfa =
       match Excitation.of_det det_i alfa with
       | Excitation.Double (_,
-            {hole = h ; particle = p ; spin = s }, 
-            {hole = h'; particle = p'; spin = s'}) -> 
+            {hole = h ; particle = p ; spin = s },
+            {hole = h'; particle = p'; spin = s'}) ->
               let fock_diag1 =
                 if s = Spin.Alfa then
                   fock_diag_alfa
@@ -1008,7 +1036,7 @@ let _pt2_en ci =
               +. (fock_diag2.(p') -. two_e h p' h p' s s'
                   +. two_e p p' p p' s s' -. two_e h' p' h' p' s' s')
       | Excitation.Single (_,
-            {hole = h ; particle = p ; spin = s }) -> 
+            {hole = h ; particle = p ; spin = s }) ->
               let fock_diag =
                 if s = Spin.Alfa then
                   fock_diag_alfa
@@ -1021,7 +1049,7 @@ let _pt2_en ci =
               h_ij mo_basis alfa alfa
       | _ -> e0.{1} -. 1.0
     in
-    let e0 = e0.{1} in 
+    let e0 = e0.{1} in
     fun alfa ->
       1. /. (e0 -. h_aa alfa)
   in
@@ -1029,16 +1057,16 @@ let _pt2_en ci =
   let mo_class = mo_class ci in
   let list_holes = List.concat
     [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
-  and list_particles = List.concat 
+  and list_particles = List.concat
     [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
   in
 
-  second_order_sum  ci  list_holes  list_particles  
+  second_order_sum  ci  list_holes  list_particles
     (is_internal ci.det_space) i_o1_alfa  i_o1_alfa  w_alfa  psi0
   |> List.fold_left (+.) 0.
 
 
-let pt2_en ci =
+let _pt2_en ci =
 
   let mo_basis  = Ds.mo_basis ci.det_space in
   let psi0, e0 = Parallel.broadcast ci.eigensystem in
@@ -1048,10 +1076,10 @@ let pt2_en ci =
   let mo_class = mo_class ci in
   let list_holes = List.concat
     [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
-  and list_particles = List.concat 
+  and list_particles = List.concat
     [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
   in
-  second_order_sum2 ci list_holes  list_particles i_o1_alfa e0.{1} psi0 
+  second_order_sum2 ci list_holes  list_particles i_o1_alfa e0.{1} psi0
 
 
 
@@ -1075,12 +1103,12 @@ let pt2_mp ci =
   let mo_class = mo_class ci in
   let list_holes = List.concat
     [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
-  and list_particles = List.concat 
+  and list_particles = List.concat
     [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
   in
 
   let psi0, _ = Parallel.broadcast ci.eigensystem in
-  second_order_sum  ci  list_holes  list_particles  
+  second_order_sum  ci  list_holes  list_particles
     (is_internal ci.det_space) i_o1_alfa  i_o1_alfa  w_alfa  psi0
   |> List.fold_left (+.) 0.
 
@@ -1097,11 +1125,11 @@ let variance ci =
   let mo_class = mo_class ci in
   let list_holes = List.concat
     [ MOClass.inactive_mos mo_class ; MOClass.active_mos mo_class ]
-  and list_particles = List.concat 
+  and list_particles = List.concat
     [ MOClass.active_mos mo_class ; MOClass.virtual_mos mo_class ]
   in
 
-  second_order_sum  ci  list_holes  list_particles  
+  second_order_sum  ci  list_holes  list_particles
     (is_internal ci.det_space) i_o1_alfa  i_o1_alfa  w_alfa  psi0
   |> List.fold_left (+.) 0.
 
@@ -1117,44 +1145,70 @@ let pt2_en_reference ci =
   let ds =
     DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis
   in
+
+  let det_stream =
+    let e =
+      let f = is_external ci.det_space  in
+      function
+      | None -> false
+      | Some d -> f d
+    in
+    let stream =
+      DeterminantSpace.determinant_stream ds
+    in
+    let rec next i =
+      let det =
+        try
+          Some (Stream.next stream)
+        with Stream.Failure -> None
+      in
+      if det = None then
+        None
+      else
+        if e det then
+          det
+        else
+          (next [@tailcall]) i
+    in
+    Stream.from next
+  in
+
   let out_dets =
-    ds
-    |> DeterminantSpace.determinants_array                                                    
-    |> Array.to_list
-    |> List.filter (fun i -> not (is_internal ci.det_space i))
+    det_stream
+    |> stream_to_list
     |> Array.of_list
   in
-      
+
   let in_dets =
     DeterminantSpace.determinants_array ci.det_space
   in
-    
+
   let m_H_aux =
     let h_aa =
       Array.map (fun ki -> h_ij aux_basis ki ki) out_dets
     in
     Array.map (fun ki ->
       Array.mapi (fun j kj ->
-        (h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j)) 
+        (h_ij aux_basis ki kj) /. (e0.{1} -. h_aa.(j))
       ) out_dets
     ) in_dets
     |> Mat.of_array
   in
-    
+
   let m_F_aux =
     Array.map (fun ki ->
       Array.map (fun kj ->
         h_ij aux_basis ki kj
-      ) out_dets 
+      ) out_dets
     ) in_dets
     |> Mat.of_array
-  in 
-    
+  in
+
   let m_HF =
     gemm m_H_aux m_F_aux ~transb:`T
   in
-  (gemm ~transa:`T psi0 @@ gemm m_HF psi0).{1,1} 
-    
+  (gemm ~transa:`T psi0 @@ gemm m_HF psi0).{1,1}
+
 
 
 
diff --git a/run_cas.ml b/run_cas.ml
index ec749f5..568214f 100644
--- a/run_cas.ml
+++ b/run_cas.ml
@@ -95,10 +95,8 @@ let () =
   if Parallel.master then
     Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);
 
-(*
   let pt2 = CI.pt2_mp ci in
   Format.fprintf ppf "CAS-MP2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);
-*)
 
 
   let pt2 = CI.pt2_en ci in