diff --git a/src/bi_ort_ints/bi_ort_ints.irp.f b/src/bi_ort_ints/bi_ort_ints.irp.f
index 6884ff38..bb894b44 100644
--- a/src/bi_ort_ints/bi_ort_ints.irp.f
+++ b/src/bi_ort_ints/bi_ort_ints.irp.f
@@ -16,108 +16,6 @@ program bi_ort_ints
 !  call test_3_e_diag
 !  call test_3_e_diag_cycle1
 ! call test_3_e
-  call routine_test_one_int
+!  call routine_test_one_int
 end
 
-subroutine routine_test_one_int
- implicit none
- integer :: p,q,r,s,ii
- integer :: i,j
- i = 3
- j = 5
- double precision :: accu
- double precision, allocatable :: vec(:)
- integer, allocatable :: iorder(:)
- allocate(vec(ao_num**4),iorder(ao_num**4))
- accu = 0.d0
- ii = 0
- do p = 1, ao_num ! 
-  do q = 1, ao_num
-   do r = 1, ao_num
-    do s = 1, ao_num
-     !<ji | ji>
-     !             
-     !                                                          j j i i
-     if(dabs(mo_l_coef(s,j) * mo_l_coef(q,i) * ao_two_e_tc_tot(s,r,q,p) * mo_r_coef(p,i) * mo_r_coef(r,j)).gt.10)then 
-     write(33,'(3(F16.10,X),4(I3,X))')mo_l_coef(s,j) * mo_l_coef(q,i)* mo_r_coef(p,i) * mo_r_coef(r,j) , ao_two_e_tc_tot(s,r,q,p), mo_l_coef(s,j) * mo_l_coef(q,i) * ao_two_e_tc_tot(s,r,q,p) * mo_r_coef(p,i) * mo_r_coef(r,j) , s,q,p,r
-     endif
-     ii += 1
-     iorder(ii) = ii 
-     vec(ii) = mo_l_coef(s,j) * mo_l_coef(q,i) * ao_two_e_tc_tot(s,r,q,p) * mo_r_coef(p,i) * mo_r_coef(r,j) 
-     accu +=  mo_l_coef(s,j) * mo_l_coef(q,i) * ao_two_e_tc_tot(s,r,q,p) * mo_r_coef(p,i) * mo_r_coef(r,j) 
-    enddo
-   enddo
-  enddo
- enddo
- call dsort(vec,iorder,ao_num**4)
- accu = 0.d0
- do i = 1, ao_num**4
-   accu += vec(i)
-   write(34,*)i,vec(i),accu
- enddo
-
- print*,'accu = ',accu
- 
-
-end
-
-subroutine routine_twoe
- implicit none
- integer :: i,j,k,l
- double precision :: old, get_mo_two_e_integral_tc_int
- double precision :: ref,new, accu, contrib, bi_ortho_mo_ints
- accu = 0.d0
- print*,'Testing the bi ortho two e'
- do j = 1, mo_num
-  do i = 1, mo_num
-   do l = 1, mo_num
-    do k = 1, mo_num
-     ! mo_non_hermit_term(k,l,i,j) = <k l| V(r_12) |i j>
-!      ref = bi_ortho_mo_ints(k,l,i,j)
-      ref = bi_ortho_mo_ints(l,k,j,i)
-      new = mo_bi_ortho_tc_two_e(l,k,j,i)
-!      old = get_mo_two_e_integral_tc_int(k,l,i,j,mo_integrals_tc_int_map)
-!      old += mo_non_hermit_term(l,k,j,i)
-
-      contrib = dabs(ref - new)
-      if(dabs(ref).gt.1.d-10)then
-       if(contrib.gt.1.d-10)then
-        print*,k,l,i,j
-        print*,ref,new,contrib
-       endif
-      endif
-      accu += contrib
-    enddo
-   enddo
-  enddo
- enddo
- print*,'accu = ',accu/(dble(mo_num)**4)
-
-end
-
-subroutine routine_onee
- implicit none
- integer :: i,k
- double precision :: ref,new,accu,contrib
- print*,'Testing the bi ortho one e'
- accu = 0.d0
- do i = 1, mo_num
-  do k = 1, mo_num
-   ref = mo_bi_ortho_tc_one_e_slow(k,i)
-   new = mo_bi_ortho_tc_one_e(k,i)
-   contrib = dabs(ref - new)
-   if(dabs(ref).gt.1.d-10)then
-    if(contrib .gt. 1.d-10)then
-     print*,'i,k',i,k
-     print*,ref,new,contrib
-    endif
-   endif
-   accu += contrib
-  enddo
- enddo
- print*,'accu = ',accu/mo_num**2
-end
-
-
-
-
diff --git a/src/fci_tc_bi/EZFIO.cfg b/src/fci_tc_bi/EZFIO.cfg
new file mode 100644
index 00000000..a2552c74
--- /dev/null
+++ b/src/fci_tc_bi/EZFIO.cfg
@@ -0,0 +1,17 @@
+[energy]
+type: double precision
+doc: Calculated Selected |FCI| energy
+interface: ezfio
+size: (determinants.n_states)
+
+[energy_pt2]
+type: double precision
+doc: Calculated |FCI| energy + |PT2|
+interface: ezfio
+size: (determinants.n_states)
+
+[cipsi_tc]
+type: character*(32)
+doc: TODO
+interface: ezfio,provider,ocaml
+default: h_tc
diff --git a/src/fci_tc_bi/NEED b/src/fci_tc_bi/NEED
new file mode 100644
index 00000000..000b0deb
--- /dev/null
+++ b/src/fci_tc_bi/NEED
@@ -0,0 +1,3 @@
+tc_bi_ortho
+davidson_undressed
+cipsi_tc_bi_ortho
diff --git a/src/fci_tc_bi/class.irp.f b/src/fci_tc_bi/class.irp.f
new file mode 100644
index 00000000..b4a68ac2
--- /dev/null
+++ b/src/fci_tc_bi/class.irp.f
@@ -0,0 +1,12 @@
+ BEGIN_PROVIDER [ logical, do_only_1h1p ]
+&BEGIN_PROVIDER [ logical, do_only_cas  ]
+&BEGIN_PROVIDER [ logical, do_ddci ]
+ implicit none
+ BEGIN_DOC
+ ! In the FCI case, all those are always false
+ END_DOC
+ do_only_1h1p = .False.
+ do_only_cas  = .False.
+ do_ddci = .False.
+END_PROVIDER
+
diff --git a/src/fci_tc_bi/copy_wf.irp.f b/src/fci_tc_bi/copy_wf.irp.f
new file mode 100644
index 00000000..cdb1ead8
--- /dev/null
+++ b/src/fci_tc_bi/copy_wf.irp.f
@@ -0,0 +1,215 @@
+
+use bitmasks
+
+subroutine copy_H_apply_buffer_to_wf_tc
+  use omp_lib
+  implicit none
+  BEGIN_DOC
+! Copies the H_apply buffer to psi_coef.
+! After calling this subroutine, N_det, psi_det and psi_coef need to be touched
+  END_DOC
+  integer(bit_kind), allocatable :: buffer_det(:,:,:)
+  double precision, allocatable  :: buffer_r_coef(:,:), buffer_l_coef(:,:)
+  integer                        :: i,j,k
+  integer                        :: N_det_old
+
+  PROVIDE H_apply_buffer_allocated
+
+
+  ASSERT (N_int > 0)
+  ASSERT (N_det > 0)
+
+  allocate ( buffer_det(N_int,2,N_det), buffer_r_coef(N_det,N_states), buffer_l_coef(N_det,N_states)  )
+
+  ! Backup determinants
+  j=0
+  do i=1,N_det
+!    if (pruned(i)) cycle  ! Pruned determinants
+    j+=1
+    ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
+    ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num)
+    buffer_det(:,:,j) = psi_det(:,:,i)
+  enddo
+  N_det_old = j
+
+  ! Backup coefficients
+  do k=1,N_states
+    j=0
+    do i=1,N_det
+!      if (pruned(i)) cycle  ! Pruned determinants
+      j += 1
+      buffer_r_coef(j,k) = psi_r_coef_bi_ortho(i,k)
+      buffer_l_coef(j,k) = psi_l_coef_bi_ortho(i,k)
+    enddo
+    ASSERT ( j == N_det_old )
+  enddo
+
+  ! Update N_det
+  N_det = N_det_old
+  do j=0,nproc-1
+    N_det = N_det + H_apply_buffer(j)%N_det
+  enddo
+
+  ! Update array sizes
+  if (psi_det_size < N_det) then
+    psi_det_size = N_det
+    TOUCH psi_det_size
+  endif
+
+  ! Restore backup in resized array
+  do i=1,N_det_old
+    psi_det(:,:,i) = buffer_det(:,:,i)
+    ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
+    ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num )
+  enddo
+  do k=1,N_states
+    do i=1,N_det_old
+      psi_r_coef_bi_ortho(i,k) = buffer_r_coef(i,k)
+      psi_l_coef_bi_ortho(i,k) = buffer_l_coef(i,k)
+    enddo
+  enddo
+
+  ! Copy new buffers
+
+  !$OMP PARALLEL DEFAULT(SHARED)                                     &
+      !$OMP PRIVATE(j,k,i) FIRSTPRIVATE(N_det_old)                   &
+      !$OMP SHARED(N_int,H_apply_buffer,psi_det,psi_r_coef_bi_ortho,psi_l_coef_bi_ortho,N_states,psi_det_size)
+  j=0
+  !$ j=omp_get_thread_num()
+  do k=0,j-1
+    N_det_old += H_apply_buffer(k)%N_det
+  enddo
+  do i=1,H_apply_buffer(j)%N_det
+    do k=1,N_int
+      psi_det(k,1,i+N_det_old) = H_apply_buffer(j)%det(k,1,i)
+      psi_det(k,2,i+N_det_old) = H_apply_buffer(j)%det(k,2,i)
+    enddo
+    ASSERT (sum(popcnt(psi_det(:,1,i+N_det_old))) == elec_alpha_num)
+    ASSERT (sum(popcnt(psi_det(:,2,i+N_det_old))) == elec_beta_num )
+  enddo
+  do k=1,N_states
+    do i=1,H_apply_buffer(j)%N_det
+      psi_r_coef_bi_ortho(i+N_det_old,k) = H_apply_buffer(j)%coef(i,k)
+      psi_l_coef_bi_ortho(i+N_det_old,k) = 0.d0
+    enddo
+  enddo
+  !$OMP BARRIER
+  H_apply_buffer(j)%N_det = 0
+  !$OMP END PARALLEL
+  SOFT_TOUCH N_det psi_det psi_r_coef_bi_ortho psi_l_coef_bi_ortho 
+
+  logical :: found_duplicates
+  call remove_duplicates_in_psi_det_tc(found_duplicates)
+  call bi_normalize(psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,size(psi_l_coef_bi_ortho,1),N_det,N_states)
+  SOFT_TOUCH N_det psi_det psi_r_coef_bi_ortho psi_l_coef_bi_ortho
+
+end
+
+subroutine remove_duplicates_in_psi_det_tc(found_duplicates)
+  implicit none
+  logical, intent(out) :: found_duplicates
+  BEGIN_DOC
+! Removes duplicate determinants in the wave function.
+  END_DOC
+  integer                        :: i,j,k
+  integer(bit_kind), allocatable :: bit_tmp(:)
+  logical,allocatable            :: duplicate(:)
+  logical                        :: dup
+
+  allocate (duplicate(N_det), bit_tmp(N_det))
+
+  found_duplicates = .False.
+
+  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,dup)
+
+  !$OMP DO
+  do i=1,N_det
+    integer, external            :: det_search_key
+    !$DIR FORCEINLINE
+    bit_tmp(i) = det_search_key(psi_det_sorted_bit_tc(1,1,i),N_int)
+    duplicate(i) = .False.
+  enddo
+  !$OMP END DO
+
+  !$OMP DO schedule(dynamic,1024)
+  do i=1,N_det-1
+    if (duplicate(i)) then
+      cycle
+    endif
+    j = i+1
+    do while (bit_tmp(j)==bit_tmp(i))
+      if (duplicate(j)) then
+        j = j+1
+        if (j > N_det) then
+          exit
+        else
+          cycle
+        endif
+      endif
+      dup = .True.
+      do k=1,N_int
+        if ( (psi_det_sorted_bit_tc(k,1,i) /= psi_det_sorted_bit_tc(k,1,j) ) &
+        .or. (psi_det_sorted_bit_tc(k,2,i) /= psi_det_sorted_bit_tc(k,2,j) ) ) then
+          dup = .False.
+          exit
+        endif
+      enddo
+      if (dup) then
+        duplicate(j) = .True.
+        found_duplicates = .True.
+      endif
+      j += 1
+      if (j > N_det) then
+        exit
+      endif
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  if (found_duplicates) then
+    k=0
+    do i=1,N_det
+      if (.not.duplicate(i)) then
+        k += 1
+        psi_det(:,:,k) = psi_det_sorted_bit_tc (:,:,i)
+        psi_r_coef_bi_ortho(k,:)  = psi_r_coef_sorted_bit(i,:)
+        psi_l_coef_bi_ortho(k,:)  = psi_l_coef_sorted_bit(i,:)
+      else
+        if (sum(abs(psi_r_coef_sorted_bit(i,:))) /= 0.d0 ) then
+          psi_r_coef_bi_ortho(k,:)  = psi_r_coef_sorted_bit(i,:)
+          psi_l_coef_bi_ortho(k,:)  = psi_l_coef_sorted_bit(i,:)
+        endif
+      endif
+    enddo
+    N_det = k
+    psi_det_sorted_bit_tc(:,:,1:N_det) = psi_det(:,:,1:N_det)
+    psi_r_coef_sorted_bit(1:N_det,:) = psi_r_coef_bi_ortho(1:N_det,:)
+    psi_l_coef_sorted_bit(1:N_det,:) = psi_l_coef_bi_ortho(1:N_det,:)
+    TOUCH N_det psi_det psi_det_sorted_bit_tc c0_weight psi_r_coef_sorted_bit psi_l_coef_sorted_bit
+  endif
+  psi_det = psi_det_sorted_tc
+  psi_r_coef_bi_ortho = psi_r_coef_sorted_bi_ortho
+  psi_l_coef_bi_ortho = psi_l_coef_sorted_bi_ortho
+  SOFT_TOUCH psi_det psi_r_coef_bi_ortho psi_l_coef_bi_ortho psi_det_sorted_bit_tc psi_r_coef_sorted_bit psi_l_coef_sorted_bit
+  deallocate (duplicate,bit_tmp)
+end
+
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_bit_tc, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_r_coef_sorted_bit, (N_det,N_states) ]
+&BEGIN_PROVIDER [ double precision, psi_l_coef_sorted_bit, (N_det,N_states) ]
+   implicit none
+   BEGIN_DOC
+   ! Determinants on which we apply $\langle i|H|psi \rangle$ for perturbation.
+   ! They are sorted by determinants interpreted as integers. Useful
+   ! to accelerate the search of a random determinant in the wave
+   ! function.
+   END_DOC
+
+   call sort_dets_by_det_search_key(N_det, psi_det, psi_r_coef_bi_ortho, size(psi_r_coef_bi_ortho,1),       &
+       psi_det_sorted_bit_tc, psi_r_coef_sorted_bit, N_states)
+   call sort_dets_by_det_search_key(N_det, psi_det, psi_l_coef_bi_ortho, size(psi_l_coef_bi_ortho,1),       &
+       psi_det_sorted_bit_tc, psi_l_coef_sorted_bit, N_states)
+
+END_PROVIDER
diff --git a/src/fci_tc_bi/diagonalize_ci.irp.f b/src/fci_tc_bi/diagonalize_ci.irp.f
new file mode 100644
index 00000000..56c561ac
--- /dev/null
+++ b/src/fci_tc_bi/diagonalize_ci.irp.f
@@ -0,0 +1,100 @@
+
+subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
+  use selection_types
+  implicit none
+  integer, intent(inout)              :: ndet      ! number of determinants from before 
+  double precision, intent(inout)     :: E_tc,norm ! E and norm from previous wave function 
+  type(pt2_type)  , intent(in)        :: pt2_data  ! PT2 from previous wave function 
+  logical, intent(in) :: print_pt2
+  BEGIN_DOC
+!  Replace the coefficients of the CI states by the coefficients of the
+!  eigenstates of the CI matrix
+  END_DOC
+  integer :: i,j
+ double precision :: pt2_tmp,pt1_norm,rpt2_tmp,abs_pt2
+ pt2_tmp = pt2_data % pt2(1)
+ abs_pt2 = pt2_data % variance(1)
+ pt1_norm = pt2_data % overlap(1,1)
+ rpt2_tmp = pt2_tmp/(1.d0 + pt1_norm)
+  print*,'*****'
+  print*,'New wave function information'
+  print*,'N_det tc               = ',N_det
+  print*,'norm_ground_left_right_bi_orth = ',norm_ground_left_right_bi_orth
+  print*,'eigval_right_tc = ',eigval_right_tc_bi_orth(1)
+  print*,'Ndet, E_tc = ',N_det,eigval_right_tc_bi_orth(1)
+  print*,'*****'
+  if(print_pt2)then
+   print*,'*****'
+   print*,'previous wave function info'
+   print*,'norm(before)      = ',norm
+   print*,'E(before)         = ',E_tc 
+   print*,'PT1 norm          = ',dsqrt(pt1_norm)
+   print*,'PT2               = ',pt2_tmp
+   print*,'rPT2              = ',rpt2_tmp
+   print*,'|PT2|             = ',abs_pt2
+   print*,'Positive PT2      = ',(pt2_tmp + abs_pt2)*0.5d0
+   print*,'Negative PT2      = ',(pt2_tmp - abs_pt2)*0.5d0
+   print*,'E(before) + PT2   = ',E_tc + pt2_tmp/norm
+   print*,'E(before) +rPT2   = ',E_tc + rpt2_tmp/norm
+   write(*,'(A28,X,I10,X,100(F16.8,X))')'Ndet,E,E+PT2,E+RPT2,|PT2|=',ndet,E_tc ,E_tc  + pt2_tmp/norm,E_tc  + rpt2_tmp/norm,abs_pt2
+   print*,'*****'
+  endif
+  E_tc  = eigval_right_tc_bi_orth(1)
+  norm  = norm_ground_left_right_bi_orth
+  ndet  = N_det
+  do j=1,N_states
+    do i=1,N_det
+      psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j)
+      psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j)
+      psi_coef(i,j) = dabs(psi_l_coef_bi_ortho(i,j) * psi_r_coef_bi_ortho(i,j))   
+    enddo
+  enddo
+  SOFT_TOUCH  eigval_left_tc_bi_orth  eigval_right_tc_bi_orth  leigvec_tc_bi_orth  reigvec_tc_bi_orth norm_ground_left_right_bi_orth psi_coef  psi_l_coef_bi_ortho psi_r_coef_bi_ortho
+
+
+
+  call save_tc_bi_ortho_wavefunction
+end
+
+subroutine print_CI_dressed(ndet, E_tc,norm,pt2_data,print_pt2)
+  use selection_types
+  implicit none
+  integer, intent(inout)              :: ndet      ! number of determinants from before 
+  double precision, intent(inout)     :: E_tc,norm ! E and norm from previous wave function 
+  type(pt2_type)  , intent(in)        :: pt2_data  ! PT2 from previous wave function 
+  logical, intent(in) :: print_pt2
+  BEGIN_DOC
+!  Replace the coefficients of the CI states by the coefficients of the
+!  eigenstates of the CI matrix
+  END_DOC
+  integer :: i,j
+  print*,'*****'
+  print*,'New wave function information'
+  print*,'N_det tc               = ',N_det
+  print*,'norm_ground_left_right_bi_orth = ',norm_ground_left_right_bi_orth
+  print*,'eigval_right_tc = ',eigval_right_tc_bi_orth(1)
+  print*,'Ndet, E_tc = ',N_det,eigval_right_tc_bi_orth(1)
+  print*,'*****'
+  if(print_pt2)then
+   print*,'*****'
+   print*,'previous wave function info'
+   print*,'norm(before)      = ',norm
+   print*,'E(before)         = ',E_tc
+   print*,'PT1 norm          = ',dsqrt(pt2_data % overlap(1,1))
+   print*,'E(before) + PT2   = ',E_tc + (pt2_data % pt2(1))/norm
+   print*,'PT2               = ',pt2_data % pt2(1)
+   print*,'Ndet, E_tc, E+PT2 = ',ndet,E_tc,E_tc + (pt2_data % pt2(1))/norm,dsqrt(pt2_data % overlap(1,1))
+   print*,'*****'
+  endif
+  E_tc  = eigval_right_tc_bi_orth(1)
+  norm  = norm_ground_left_right_bi_orth
+  ndet  = N_det
+  do j=1,N_states
+    do i=1,N_det
+      psi_coef(i,j) = reigvec_tc_bi_orth(i,j)
+    enddo
+  enddo
+  SOFT_TOUCH  eigval_left_tc_bi_orth  eigval_right_tc_bi_orth  leigvec_tc_bi_orth  norm_ground_left_right_bi_orth  psi_coef  reigvec_tc_bi_orth 
+
+end
+
diff --git a/src/fci_tc_bi/fci_tc_bi_ortho.irp.f b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f
new file mode 100644
index 00000000..84ac8166
--- /dev/null
+++ b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f
@@ -0,0 +1,85 @@
+program fci
+  implicit none
+  BEGIN_DOC
+  ! Selected Full Configuration Interaction with stochastic selection
+  ! and PT2.
+  !
+  ! This program performs a |CIPSI|-like selected |CI| using a
+  ! stochastic scheme for both the selection of the important Slater
+  ! determinants and the computation of the |PT2| correction. This
+  ! |CIPSI|-like algorithm will be performed for the lowest states of
+  ! the variational space (see :option:`determinants n_states`). The
+  ! |FCI| program will stop when reaching at least one the two following
+  ! conditions:
+  !
+  ! * number of Slater determinants > :option:`determinants n_det_max`
+  ! * abs(|PT2|) less than :option:`perturbation pt2_max`
+  !
+  ! The following other options can be of interest:
+  !
+  ! :option:`determinants read_wf`
+  !   When set to |false|, the program starts with a ROHF-like Slater
+  !   determinant as a guess wave function. When set to |true|, the
+  !   program starts with the wave function(s) stored in the |EZFIO|
+  !   directory as guess wave function(s).
+  !
+  ! :option:`determinants s2_eig`
+  !   When set to |true|, the selection will systematically add all the
+  !   necessary Slater determinants in order to have a pure spin wave
+  !   function with an |S^2| value corresponding to
+  !   :option:`determinants expected_s2`.
+  !
+  ! For excited states calculations, it is recommended to start with
+  ! :ref:`cis` or :ref:`cisd` guess wave functions, eventually in
+  ! a restricted set of |MOs|, and to set :option:`determinants s2_eig`
+  ! to |true|.
+  !
+  END_DOC
+
+
+  my_grid_becke = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid 
+  pruning = -1.d0
+  touch pruning
+!  pt2_relative_error = 0.01d0
+!  touch pt2_relative_error
+  call run_cipsi_tc
+
+end
+
+
+subroutine run_cipsi_tc
+
+  implicit none
+
+  if (.not.is_zmq_slave) then
+    PROVIDE psi_det psi_coef mo_bi_ortho_tc_two_e mo_bi_ortho_tc_one_e
+    if(elec_alpha_num+elec_beta_num.ge.3)then
+      if(three_body_h_tc)then
+        call provide_all_three_ints_bi_ortho
+      endif
+    endif
+    ! ---
+
+    if (do_pt2) then
+      call run_stochastic_cipsi
+    else
+      call run_cipsi
+    endif
+
+  else
+    PROVIDE mo_bi_ortho_tc_one_e mo_bi_ortho_tc_two_e pt2_min_parallel_tasks
+    if(elec_alpha_num+elec_beta_num.ge.3)then
+      if(three_body_h_tc)then
+        call provide_all_three_ints_bi_ortho
+      endif
+    endif
+    ! ---
+
+    call run_slave_cipsi
+
+  endif
+
+end
diff --git a/src/fci_tc_bi/generators.irp.f b/src/fci_tc_bi/generators.irp.f
new file mode 100644
index 00000000..55c0cbb9
--- /dev/null
+++ b/src/fci_tc_bi/generators.irp.f
@@ -0,0 +1,51 @@
+use bitmasks
+
+BEGIN_PROVIDER [ integer, N_det_generators ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the number of generators is 1 : the
+ ! Hartree-Fock determinant
+ END_DOC
+ integer :: i
+ double precision :: norm
+ call write_time(6)
+ norm = 1.d0
+ N_det_generators = N_det
+ do i=1,N_det
+   norm = norm - psi_average_norm_contrib_sorted_tc(i)
+   if (norm - 1.d-10 < 1.d0 - threshold_generators) then
+     N_det_generators = i
+     exit
+   endif
+ enddo
+ N_det_generators = max(N_det_generators,1)
+ call write_int(6,N_det_generators,'Number of generators')
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted_tc(1:N_int,1:2,1:N_det)
+ psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted_tc(1:N_det,1:N_states)
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_gen, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc_gen, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer, psi_det_sorted_tc_gen_order,     (psi_det_size)  ]
+
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_det_sorted_tc_gen = psi_det_sorted_tc
+ psi_coef_sorted_tc_gen = psi_coef_sorted_tc
+ psi_det_sorted_tc_gen_order = psi_det_sorted_tc_order
+END_PROVIDER
+
+
diff --git a/src/fci_tc_bi/save_energy.irp.f b/src/fci_tc_bi/save_energy.irp.f
new file mode 100644
index 00000000..7c41d00f
--- /dev/null
+++ b/src/fci_tc_bi/save_energy.irp.f
@@ -0,0 +1,9 @@
+subroutine save_energy(E,pt2)
+  implicit none
+  BEGIN_DOC
+! Saves the energy in |EZFIO|.
+  END_DOC
+  double precision, intent(in) :: E(N_states), pt2(N_states)
+  call ezfio_set_fci_tc_energy(E(1:N_states))
+  call ezfio_set_fci_tc_energy_pt2(E(1:N_states)+pt2(1:N_states))
+end
diff --git a/src/fci_tc_bi/scripts_fci_tc/CH2.xyz b/src/fci_tc_bi/scripts_fci_tc/CH2.xyz
new file mode 100644
index 00000000..9fa57f4b
--- /dev/null
+++ b/src/fci_tc_bi/scripts_fci_tc/CH2.xyz
@@ -0,0 +1,6 @@
+3
+
+C                    6.000000     0.000000     0.000000     0.173480    
+H                    1.000000     0.000000    -0.861500    -0.520430
+H                    1.000000     0.000000     0.861500    -0.520430
+
diff --git a/src/fci_tc_bi/scripts_fci_tc/FH.xyz b/src/fci_tc_bi/scripts_fci_tc/FH.xyz
new file mode 100644
index 00000000..9a1563f4
--- /dev/null
+++ b/src/fci_tc_bi/scripts_fci_tc/FH.xyz
@@ -0,0 +1,5 @@
+2 
+
+H                    0.000000     0.000000    -0.825120
+F                    0.000000     0.000000     0.091680
+
diff --git a/src/fci_tc_bi/scripts_fci_tc/extract_tables.sh b/src/fci_tc_bi/scripts_fci_tc/extract_tables.sh
new file mode 100755
index 00000000..a585884e
--- /dev/null
+++ b/src/fci_tc_bi/scripts_fci_tc/extract_tables.sh
@@ -0,0 +1,16 @@
+
+input=h2o
+basis=dz
+EZFIO=${input}_${basis}_bi_ortho
+file=${EZFIO}.tc_fci.out 
+grep "Ndet,E,E+PT2,E+RPT2,|PT2|=" ${file} | cut -d "=" -f 2 > data_${EZFIO}
+file=${EZFIO}.tc_fci_normal_order.out
+grep "Ndet,E,E+PT2,E+RPT2=" ${file} | cut -d "=" -f 2 > data_${EZFIO}_normal
+
+#EZFIO=${input}_${basis}_ortho
+#file=${EZFIO}.tc_fci.out 
+#grep "Ndet, E_tc, E+PT2 =" ${file} | cut -d "=" -f 2 > data_${EZFIO}
+#file=${EZFIO}.tc_fci_normal_order.out
+#grep "Ndet, E_tc, E+PT2 =" ${file} | cut -d "=" -f 2 > data_${EZFIO}_normal
+
+#zip data_${input}_${basis}.zip data*
diff --git a/src/fci_tc_bi/scripts_fci_tc/h2o.sh b/src/fci_tc_bi/scripts_fci_tc/h2o.sh
new file mode 100644
index 00000000..d0afca30
--- /dev/null
+++ b/src/fci_tc_bi/scripts_fci_tc/h2o.sh
@@ -0,0 +1,41 @@
+#!/bin/bash
+# This is a sample PBS script
+# temps CPU a ajuster au calcul
+   #PBS -l cput=2000:00:00
+   #PBS -l nodes=1:ppn=16
+# memoire a ajuster au calcul
+   #PBS -l vmem=100gb
+# a changer
+# Pour savoir sur quel noeud on est
+#echo $HOSTNAME
+# Startdir = ou sont les fichiers d'input, par defaut HOMEDIR
+#
+StartDir=$PBS_O_WORKDIR
+echo $StartDir
+#
+# SCRATCHDIR = espace temporaire (local au noeud et a vider apres le calcul)
+# NE PAS MODIFIER
+ulimit -s unlimited
+export SCRATCHDIR=/scratch/$USER/$PBS_JOBID
+#
+cd $StartDir
+
+
+############################################################################
+#### EXAMPLE OF SCRIPT TO RUN A CIPSI CALCULATION ON 5 STATES ON THE Ne^+ CATION
+#### USING NATURAL ORBITALS OF A SMALL CIPSI AS MOS 
+#### ALL STATES WILL HAVE THE SAME SPIN SIMETRY : A DOUBLET 
+
+####### YOU PUT THE PATH TO YOUR 
+QP_ROOT=/home_lct/eginer/programs/qp2
+source ${QP_ROOT}/quantum_package.rc 
+####### YOU LOAD SOME LIBRARIES 
+alias python3='/programmes/installation/Python/3.7.1/bin/python3'
+type -a python3
+
+export OMP_NUM_THREADS=16
+
+module load intel2016_OMPI-V2
+
+source ~/programs/qp2/quantum_package.rc
+./script.sh h2o dz O 1
diff --git a/src/fci_tc_bi/scripts_fci_tc/h2o.xyz b/src/fci_tc_bi/scripts_fci_tc/h2o.xyz
new file mode 100644
index 00000000..dee51ffc
--- /dev/null
+++ b/src/fci_tc_bi/scripts_fci_tc/h2o.xyz
@@ -0,0 +1,6 @@
+3
+
+O      0.000000     0.000000     0.000000
+H      0.000000     0.000000     0.957200
+H     -0.926627     0.000000    -0.239987
+
diff --git a/src/fci_tc_bi/scripts_fci_tc/script.sh b/src/fci_tc_bi/scripts_fci_tc/script.sh
new file mode 100755
index 00000000..58585658
--- /dev/null
+++ b/src/fci_tc_bi/scripts_fci_tc/script.sh
@@ -0,0 +1,31 @@
+source /home_lct/eginer/qp2/quantum_package.rc
+input=$1
+ basis=$2
+ atom=$3
+ mul=$4
+ EXPORT_OMP_NUM_THREADS=16
+ dir=${input}_${basis}
+ mkdir ${dir}
+ cp ${input}.xyz ${dir}/
+ cd $dir
+ EZFIO=${input}_${basis}_bi_ortho
+ qp create_ezfio -b "${atom}:cc-pcv${basis}|H:cc-pv${basis}" ${input}.xyz -m $mul -o $EZFIO
+ qp run scf
+ # Getting THE GOOD VALUE OF MU
+ qp run print_mu_av_tc  | tee ${EZFIO_FILE}.mu_av.out
+ mu=`grep "average_mu_rs_c_lda  =" ${EZFIO_FILE}.mu_av.out | cut -d "=" -f 2`
+ qp set ao_two_e_erf_ints mu_erf $mu
+ # Carrying the BI-ORTHO TC-SCF 
+ qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
+ #   Three body terms without normal order
+ ### THREE E TERMS FOR FCI
+ qp set tc_keywords three_body_h_tc True  
+ qp set tc_keywords double_normal_ord False
+ qp set perturbation pt2_max 0.003
+ qp run fci_tc_bi_ortho | tee ${EZFIO_FILE}.tc_fci.out
+ #   Three body terms with normal order
+ qp set tc_keywords double_normal_ord True
+ qp run fci_tc_bi_ortho | tee ${EZFIO_FILE}.tc_fci_normal_order.out
+
+cd ../
+
diff --git a/src/fci_tc_bi/selectors.irp.f b/src/fci_tc_bi/selectors.irp.f
new file mode 100644
index 00000000..734c8ed0
--- /dev/null
+++ b/src/fci_tc_bi/selectors.irp.f
@@ -0,0 +1,92 @@
+use bitmasks
+
+BEGIN_PROVIDER [ double precision, threshold_selectors ]
+ implicit none
+ BEGIN_DOC
+ ! Thresholds on selectors (fraction of the square of the norm)
+ END_DOC
+ threshold_selectors = dsqrt(threshold_generators)
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, N_det_selectors]
+  implicit none
+  BEGIN_DOC
+  ! For Single reference wave functions, the number of selectors is 1 : the
+  ! Hartree-Fock determinant
+  END_DOC
+  integer                        :: i
+  double precision               :: norm, norm_max
+  call write_time(6)
+  N_det_selectors = N_det
+  norm = 1.d0
+  do i=1,N_det
+    norm = norm - psi_average_norm_contrib_tc(i)
+    if (norm - 1.d-10 < 1.d0 - threshold_selectors) then
+      N_det_selectors = i
+      exit
+    endif
+  enddo
+  N_det_selectors = max(N_det_selectors,N_det_generators)
+  call write_int(6,N_det_selectors,'Number of selectors')
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_selectors, (N_int,2,psi_selectors_size) ]
+&BEGIN_PROVIDER [ double precision, psi_selectors_coef, (psi_selectors_size,N_states) ]
+  implicit none
+  BEGIN_DOC
+  ! Determinants on which we apply <i|H|psi> for perturbation.
+  END_DOC
+  integer                        :: i,k
+
+  do i=1,N_det_selectors
+    do k=1,N_int
+      psi_selectors(k,1,i) = psi_det_sorted_tc(k,1,i)
+      psi_selectors(k,2,i) = psi_det_sorted_tc(k,2,i)
+    enddo
+  enddo
+  do k=1,N_states
+    do i=1,N_det_selectors
+      psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k)
+    enddo
+  enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp, (N_states,psi_selectors_size) ]
+  implicit none
+  BEGIN_DOC
+  ! Transposed psi_selectors
+  END_DOC
+  integer                        :: i,k
+
+  do i=1,N_det_selectors
+    do k=1,N_states
+      psi_selectors_coef_transp(k,i) = psi_selectors_coef(i,k)
+    enddo
+  enddo
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, psi_selectors_rcoef_bi_orth_transp, (N_states, psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_selectors_lcoef_bi_orth_transp, (N_states, psi_det_size) ]
+
+  implicit none
+  integer :: i, k
+
+  psi_selectors_rcoef_bi_orth_transp = 0.d0
+  psi_selectors_lcoef_bi_orth_transp = 0.d0
+
+  print*,'N_det,N_det_selectors',N_det,N_det_selectors
+  do i = 1, N_det_selectors
+    do k = 1, N_states
+      psi_selectors_rcoef_bi_orth_transp(k,i) = psi_r_coef_sorted_bi_ortho(i,k)
+      psi_selectors_lcoef_bi_orth_transp(k,i) = psi_l_coef_sorted_bi_ortho(i,k)
+    enddo
+  enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, psi_selectors_size ]
+ implicit none
+ psi_selectors_size = psi_det_size
+END_PROVIDER
+
diff --git a/src/fci_tc_bi/zmq.irp.f b/src/fci_tc_bi/zmq.irp.f
new file mode 100644
index 00000000..cb2df483
--- /dev/null
+++ b/src/fci_tc_bi/zmq.irp.f
@@ -0,0 +1,103 @@
+BEGIN_TEMPLATE
+
+integer function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Put $X on the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer                        :: rc
+  character*(256)                :: msg
+
+  zmq_put_$X = 0
+
+  write(msg,'(A,1X,I8,1X,A200)') 'put_data '//trim(zmq_state), worker_id, '$X'
+  rc = f77_zmq_send(zmq_to_qp_run_socket,trim(msg),len(trim(msg)),ZMQ_SNDMORE)
+  if (rc /= len(trim(msg))) then
+    zmq_put_$X = -1
+    return
+  endif
+
+  rc = f77_zmq_send(zmq_to_qp_run_socket,$X,4,0)
+  if (rc /= 4) then
+    zmq_put_$X = -1
+    return
+  endif
+
+  rc = f77_zmq_recv(zmq_to_qp_run_socket,msg,len(msg),0)
+  if (msg(1:rc) /= 'put_data_reply ok') then
+    zmq_put_$X = -1
+    return
+  endif
+
+end
+
+integer function zmq_get_$X(zmq_to_qp_run_socket, worker_id)
+  use f77_zmq
+  implicit none
+  BEGIN_DOC
+! Get $X from the qp_run scheduler
+  END_DOC
+  integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
+  integer, intent(in)            :: worker_id
+  integer                        :: rc
+  character*(256)                :: msg
+
+  PROVIDE zmq_state
+  zmq_get_$X = 0
+  if (mpi_master) then
+
+    write(msg,'(A,1X,I8,1X,A200)') 'get_data '//trim(zmq_state), worker_id, '$X'
+    rc = f77_zmq_send(zmq_to_qp_run_socket,trim(msg),len(trim(msg)),0)
+    if (rc /= len(trim(msg))) then
+      zmq_get_$X = -1
+      go to 10
+    endif
+
+    rc = f77_zmq_recv(zmq_to_qp_run_socket,msg,len(msg),0)
+    if (msg(1:14) /= 'get_data_reply') then
+      zmq_get_$X = -1
+      go to 10
+    endif
+
+    rc = f77_zmq_recv(zmq_to_qp_run_socket,$X,4,0)
+    if (rc /= 4) then
+      zmq_get_$X = -1
+      go to 10
+    endif
+
+  endif
+
+  10 continue
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+    include 'mpif.h'
+    integer :: ierr
+
+    call MPI_BCAST (zmq_get_$X, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      print *,  irp_here//': Unable to broadcast N_det_generators'
+      stop -1
+    endif
+    call MPI_BCAST ($X, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      print *,  irp_here//': Unable to broadcast N_det_generators'
+      stop -1
+    endif
+  IRP_ENDIF
+
+end
+
+SUBST [ X ]
+
+N_det_generators ;;
+N_det_selectors ;;
+
+END_TEMPLATE
+