diff --git a/src/Davidson/davidson_slave.irp.f b/src/Davidson/davidson_slave.irp.f
new file mode 100644
index 00000000..d8143958
--- /dev/null
+++ b/src/Davidson/davidson_slave.irp.f
@@ -0,0 +1,32 @@
+program davidson_slave
+  use f77_zmq
+  implicit none
+
+
+  integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
+  integer(ZMQ_PTR) :: zmq_to_qp_run_socket
+  double precision :: energy(N_states_diag)
+  character*(64) :: state
+  
+  call provide_everything
+  call switch_qp_run_to_master
+  call omp_set_nested(.True.)
+  
+  zmq_context = f77_zmq_ctx_new ()
+  zmq_state = 'davidson'
+  state = 'Waiting'
+
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
+  do
+    call wait_for_state(zmq_state,state)
+    if(trim(state) /= "davidson") exit
+    integer :: rc, i
+    print *,  'Davidson slave running'
+    call davidson_slave_tcp(i)
+  end do
+end
+
+subroutine provide_everything
+  PROVIDE mo_bielec_integrals_in_map psi_det_sorted_bit N_states_diag zmq_context ref_bitmask_energy
+end subroutine
+
diff --git a/src/Davidson/diag_restart_save_all_states.irp.f b/src/Davidson/diag_restart_save_all_states.irp.f
new file mode 100644
index 00000000..3bdc37c5
--- /dev/null
+++ b/src/Davidson/diag_restart_save_all_states.irp.f
@@ -0,0 +1,16 @@
+program diag_and_save
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call routine
+end
+
+subroutine routine
+ implicit none
+ call diagonalize_CI
+ print*,'N_det = ',N_det
+ call save_wavefunction_general(N_det,N_states_diag,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
+
+
+
+end
diff --git a/src/Davidson/diag_restart_save_lowest_state.irp.f b/src/Davidson/diag_restart_save_lowest_state.irp.f
new file mode 100644
index 00000000..0e379aae
--- /dev/null
+++ b/src/Davidson/diag_restart_save_lowest_state.irp.f
@@ -0,0 +1,29 @@
+program diag_and_save
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call routine
+end
+
+subroutine routine
+ implicit none
+ print*,'N_det = ',N_det
+ PROVIDE H_apply_buffer_allocated
+ if (s2_eig) then
+   call make_s2_eigenfunction
+ endif
+ call diagonalize_CI
+ integer :: igood_state
+ igood_state=1
+ double precision, allocatable :: psi_coef_tmp(:)
+ allocate(psi_coef_tmp(n_det))
+ integer :: i
+ do i = 1, N_det
+  psi_coef_tmp(i) = psi_coef(i,igood_state)
+ enddo
+ call save_wavefunction_general(N_det,1,psi_det,n_det,psi_coef_tmp)
+ deallocate(psi_coef_tmp)
+
+
+
+end
diff --git a/src/Davidson/diag_restart_save_one_state.irp.f b/src/Davidson/diag_restart_save_one_state.irp.f
new file mode 100644
index 00000000..c5f4e59d
--- /dev/null
+++ b/src/Davidson/diag_restart_save_one_state.irp.f
@@ -0,0 +1,26 @@
+program diag_and_save
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call routine
+end
+
+subroutine routine
+ implicit none
+ print*,'N_det = ',N_det
+ call diagonalize_CI
+ write(*,*)'Which state would you like to save ?'
+ integer :: igood_state
+ read(5,*)igood_state
+ double precision, allocatable :: psi_coef_tmp(:)
+ allocate(psi_coef_tmp(n_det))
+ integer :: i
+ do i = 1, N_det
+  psi_coef_tmp(i) = psi_coef(i,igood_state)
+ enddo
+ call save_wavefunction_general(N_det,1,psi_det,n_det,psi_coef_tmp)
+ deallocate(psi_coef_tmp)
+
+
+
+end
diff --git a/src/Davidson/guess_lowest_state.irp.f b/src/Davidson/guess_lowest_state.irp.f
new file mode 100644
index 00000000..f6d0a004
--- /dev/null
+++ b/src/Davidson/guess_lowest_state.irp.f
@@ -0,0 +1,162 @@
+program first_guess
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Select all the determinants with the lowest energy as a starting point.
+  END_DOC
+  integer                        :: i,j
+  double precision, allocatable  :: orb_energy(:)
+  double precision               :: E
+  integer, allocatable           :: kept(:)
+  integer                        :: nelec_kept(2)
+  character                      :: occ_char, keep_char
+  
+  PROVIDE H_apply_buffer_allocated psi_det
+  allocate (orb_energy(mo_tot_num), kept(0:mo_tot_num))
+  nelec_kept(1:2) = 0
+  kept(0) = 0
+
+  print *,  'Orbital energies'
+  print *,  '================'
+  print *,  ''
+  do i=1,mo_tot_num
+    keep_char = ' '
+    occ_char = '-'
+    orb_energy(i) = mo_mono_elec_integral(i,i)
+    do j=1,elec_beta_num
+      if (i==j) cycle
+      orb_energy(i) += mo_bielec_integral_jj_anti(i,j)
+    enddo
+    do j=1,elec_alpha_num
+      orb_energy(i) += mo_bielec_integral_jj(i,j)
+    enddo
+    if ( (orb_energy(i) > -.5d0).and.(orb_energy(i) < .1d0) ) then
+      kept(0) += 1
+      keep_char = 'X'
+      kept( kept(0) ) = i
+      if (i <= elec_beta_num) then
+        nelec_kept(2) += 1
+      endif
+      if (i <= elec_alpha_num) then
+        nelec_kept(1) += 1
+      endif
+    endif
+    if (i <= elec_alpha_num) then
+      if (i <= elec_beta_num) then
+        occ_char = '#'
+      else
+        occ_char = '+'
+      endif
+    endif
+    print '(I4, 3X, A, 3X, F10.6, 3X, A)', i, occ_char, orb_energy(i), keep_char 
+  enddo
+
+  
+  integer, allocatable           :: list  (:,:)
+  integer(bit_kind), allocatable :: string(:,:)
+  allocate ( list(N_int*bit_kind_size,2), string(N_int,2) )
+
+  string = ref_bitmask
+  call bitstring_to_list( string(1,1), list(1,1), elec_alpha_num, N_int)
+  call bitstring_to_list( string(1,2), list(1,2), elec_beta_num , N_int)
+
+  psi_det_alpha_unique(:,1) = string(:,1)
+  psi_det_beta_unique (:,1) = string(:,2)
+  N_det_alpha_unique = 1
+  N_det_beta_unique  = 1
+
+  integer :: i1,i2,i3,i4,i5,i6,i7,i8,i9
+
+  psi_det_size = kept(0)**(nelec_kept(1)+nelec_kept(2))
+  print *,  kept(0), nelec_kept(:)
+  call write_int(6,psi_det_size,'psi_det_size')
+  TOUCH psi_det_size
+
+BEGIN_SHELL [ /usr/bin/python ]
+
+template_alpha_ext = """
+do %(i2)s = %(i1)s-1,1,-1
+  list(elec_alpha_num-%(i)d,1) = kept(%(i2)s)
+  call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
+"""
+
+template_alpha = """
+do %(i2)s = %(i1)s-1,1,-1
+  list(elec_alpha_num-%(i)d,1) = kept(%(i2)s)
+  call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
+  N_det_alpha_unique += 1
+  psi_det_alpha_unique(:,N_det_alpha_unique) = string(:,1)
+"""
+
+template_beta_ext = """
+do %(i2)s = %(i1)s-1,1,-1
+  list(elec_beta_num-%(i)d,2) = kept(%(i2)s)
+  call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
+"""
+template_beta = """
+do %(i2)s = %(i1)s-1,1,-1
+  list(elec_beta_num-%(i)d,2) = kept(%(i2)s)
+  call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
+  N_det_beta_unique += 1
+  psi_det_beta_unique(:,N_det_beta_unique) = string(:,2)
+"""
+
+def write(template_ext,template,imax):
+    print "case(%d)"%(imax)
+    def aux(i2,i1,i,j):
+        if (i==imax-1):
+          print template%locals()
+        else:
+          print template_ext%locals()
+        i += 1
+        j -= 1
+        if (i != imax):
+            i1 = "i%d"%(i)
+            i2 = "i%d"%(i+1)
+            aux(i2,i1,i,j)
+        print "enddo"
+
+    i2 = "i1"
+    i1 = "kept(0)+1"
+    i  = 0
+    aux(i2,i1,i,imax)
+
+def main():
+    print """
+    select case (nelec_kept(1))
+      case(0)
+        continue
+    """
+    for imax in range(1,10):
+        write(template_alpha_ext,template_alpha,imax)
+
+    print """
+    end select
+
+    select case (nelec_kept(2))
+      case(0)
+        continue
+    """
+    for imax in range(1,10):
+        write(template_beta_ext,template_beta,imax)
+    print "end select"
+
+main()
+
+END_SHELL
+
+  TOUCH N_det_alpha_unique N_det_beta_unique psi_det_alpha_unique psi_det_beta_unique
+  call create_wf_of_psi_bilinear_matrix(.False.)
+  call diagonalize_ci
+  j= N_det
+  do i=1,N_det
+    if (psi_average_norm_contrib_sorted(i) < 1.d-6) then
+       j = i-1
+       exit
+    endif
+!  call debug_det(psi_det_sorted(1,1,i),N_int)
+  enddo
+  call save_wavefunction_general(j,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
+
+  deallocate(orb_energy, kept, list, string)
+end
diff --git a/src/Davidson/null_dressing_vector.irp.f b/src/Davidson/null_dressing_vector.irp.f
new file mode 100644
index 00000000..faffe964
--- /dev/null
+++ b/src/Davidson/null_dressing_vector.irp.f
@@ -0,0 +1,10 @@
+ BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ]
+&BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Null dressing vectors
+ END_DOC
+ dressing_column_h(:,:) = 0.d0
+ dressing_column_s(:,:) = 0.d0
+END_PROVIDER
+
diff --git a/src/Davidson/print_H_matrix_restart.irp.f b/src/Davidson/print_H_matrix_restart.irp.f
new file mode 100644
index 00000000..57fc3633
--- /dev/null
+++ b/src/Davidson/print_H_matrix_restart.irp.f
@@ -0,0 +1,176 @@
+program print_H_matrix_restart
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call routine
+
+end
+
+subroutine routine
+ use bitmasks
+ implicit none
+ integer :: i,j
+ integer, allocatable :: H_matrix_degree(:,:)
+ double precision, allocatable :: H_matrix_phase(:,:)
+ integer :: degree
+ integer(bit_kind), allocatable :: keys_tmp(:,:,:)
+ allocate(keys_tmp(N_int,2,N_det))
+ do i = 1, N_det
+  print*,''
+  call debug_det(psi_det(1,1,i),N_int)
+  do j = 1, N_int
+   keys_tmp(j,1,i) = psi_det(j,1,i)
+   keys_tmp(j,2,i) = psi_det(j,2,i)
+  enddo
+ enddo
+ if(N_det.ge.10000)then
+  print*,'Warning !!!'
+  print*,'Number of determinants is ',N_det
+  print*,'It means that the H matrix will be enormous !'
+  print*,'stoppping ..'
+  stop
+ endif
+ print*,''
+ print*,'Determinants '
+ do i = 1, N_det
+ enddo
+ allocate(H_matrix_degree(N_det,N_det),H_matrix_phase(N_det,N_det))
+ integer         :: exc(0:2,2,2)
+ double precision  :: phase
+ do i = 1, N_det
+  do j = i, N_det 
+   call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+   H_matrix_degree(i,j) = degree
+   H_matrix_degree(j,i) = degree
+   phase = 0.d0
+   if(degree==1.or.degree==2)then
+    call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
+   endif
+   H_matrix_phase(i,j) = phase
+   H_matrix_phase(j,i) = phase
+  enddo
+ enddo
+ print*,'H matrix '
+ double precision :: ref_h_matrix,s2
+ ref_h_matrix = H_matrix_all_dets(1,1)
+ print*,'HF like determinant energy = ',ref_bitmask_energy+nuclear_repulsion
+ print*,'Ref element of H_matrix    = ',ref_h_matrix+nuclear_repulsion
+ print*,'Printing the H matrix ...'
+ print*,''
+ print*,''
+!do i = 1, N_det
+! H_matrix_all_dets(i,i) -= ref_h_matrix
+!enddo
+
+ do i = 1, N_det
+  H_matrix_all_dets(i,i) += nuclear_repulsion
+ enddo
+
+!do i = 5,N_det
+! H_matrix_all_dets(i,3) = 0.d0
+! H_matrix_all_dets(3,i) = 0.d0
+! H_matrix_all_dets(i,4) = 0.d0
+! H_matrix_all_dets(4,i) = 0.d0
+!enddo
+
+
+
+
+ 
+ do i = 1, N_det
+  write(*,'(I3,X,A3,1000(F16.7))')i,' | ',H_matrix_all_dets(i,:)
+ enddo
+
+ print*,''
+ print*,''
+ print*,''
+ print*,'Printing the degree of excitations within the H matrix'
+ print*,''
+ print*,''
+ do i = 1, N_det
+  write(*,'(I3,X,A3,X,1000(I1,X))')i,' | ',H_matrix_degree(i,:)
+ enddo
+
+
+ print*,''
+ print*,''
+ print*,'Printing the phase of the Hamiltonian matrix elements '
+ print*,''
+ print*,''
+ do i = 1, N_det
+  write(*,'(I3,X,A3,X,1000(F3.0,X))')i,' | ',H_matrix_phase(i,:)
+ enddo
+ print*,''
+
+
+ double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:)
+ double precision, allocatable  :: s2_eigvalues(:)
+ allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
+ allocate (eigenvalues(N_det),s2_eigvalues(N_det))
+ call lapack_diag(eigenvalues,eigenvectors,                       &
+     H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
+ print*,'Two first eigenvectors '
+ call u_0_S2_u_0(s2_eigvalues,eigenvectors,n_det,keys_tmp,N_int,N_det,size(eigenvectors,1))
+ do j =1, N_states
+   print*,'s2 = ',s2_eigvalues(j)
+   print*,'e  = ',eigenvalues(j)
+   print*,'coefs : '
+   do i = 1, N_det
+    print*,'i = ',i,eigenvectors(i,j)
+   enddo
+   if(j>1)then
+    print*,'Delta E(H)  = ',eigenvalues(1) - eigenvalues(j)
+    print*,'Delta E(eV) = ',(eigenvalues(1) - eigenvalues(j))*27.2114d0
+   endif
+ enddo
+ double precision               :: get_mo_bielec_integral,k_a_iv,k_b_iv
+ integer :: h1,p1,h2,p2
+ h1 = 10
+ p1 = 16 
+ h2 = 14 
+ p2 = 14
+!h1 = 1
+!p1 = 4
+!h2 = 2
+!p2 = 2
+ k_a_iv = get_mo_bielec_integral(h1,h2,p2,p1,mo_integrals_map)
+ h2 = 15 
+ p2 = 15
+ k_b_iv = get_mo_bielec_integral(h1,h2,p2,p1,mo_integrals_map)
+ print*,'k_a_iv = ',k_a_iv
+ print*,'k_b_iv = ',k_b_iv
+ double precision :: k_av,k_bv,k_ai,k_bi
+ h1 = 16
+ p1 = 14 
+ h2 = 14 
+ p2 = 16
+ k_av = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
+ h1 = 16
+ p1 = 15 
+ h2 = 15 
+ p2 = 16
+ k_bv = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
+ 
+ h1 = 10
+ p1 = 14 
+ h2 = 14 
+ p2 = 10
+ k_ai = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
+ 
+ h1 = 10
+ p1 = 15 
+ h2 = 15 
+ p2 = 10
+ k_bi = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
+ 
+ print*,'k_av, k_bv = ',k_av,k_bv
+ print*,'k_ai, k_bi = ',k_ai,k_bi
+ double precision :: k_iv
+
+ h1 = 10
+ p1 = 16 
+ h2 = 16 
+ p2 = 10
+ k_iv = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
+ print*,'k_iv       = ',k_iv
+end
diff --git a/src/Davidson/print_energy.irp.f b/src/Davidson/print_energy.irp.f
new file mode 100644
index 00000000..ae6f1da2
--- /dev/null
+++ b/src/Davidson/print_energy.irp.f
@@ -0,0 +1,22 @@
+program print_energy
+ implicit none
+ read_wf = .true.
+ touch read_wf
+ call routine
+end
+
+subroutine routine
+ implicit none
+ integer :: i,j
+ double precision :: accu,hij
+
+ print*, 'psi_energy          = ',psi_energy + nuclear_repulsion
+ accu = 0.d0
+! do i = 1,N_det
+!  do j = 1,N_det
+!   call i_H_j(psi_det(1,1,j),psi_det(1,1,i),N_int,hij)
+!   accu += psi_coef(i,1) * psi_coef(j,1) * hij
+!  enddo
+! enddo
+! print*, 'accu                = ',accu + nuclear_repulsion
+end