diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml
index 17d31cea..5a656d2d 100644
--- a/ocaml/qp_run.ml
+++ b/ocaml/qp_run.ml
@@ -134,7 +134,8 @@ let run slave exe ezfio_file =
   let duration = Time.diff (Time.now()) time_start 
   |> Core.Span.to_string in
   Printf.printf "Wall time : %s\n\n" duration;
-  exit exit_code
+  if (exit_code <> 0) then
+    exit exit_code
 
 let spec = 
   let open Command.Spec in
diff --git a/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg
index 8e4b2847..df55c554 100644
--- a/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg
+++ b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg
@@ -10,5 +10,17 @@ doc: Exponents of the additional Slater functions
 size: (nuclei.nucl_num,mo_basis.mo_tot_num)
 interface: ezfio, provider
 
+[projector]
+type: double precision
+doc: Orthogonal AO basis
+size: (ao_basis.ao_num,ao_basis.ao_num)
+interface: ezfio
+
+[ao_orthoSlaOverlap]
+type: double precision
+doc: Orthogonal AO basis
+size: (ao_basis.ao_num,nuclei.nucl_num)
+interface: ezfio
+
 
 
diff --git a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f
index 490e344c..61a5a49c 100644
--- a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f
+++ b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f
@@ -13,7 +13,7 @@ BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ]
       cusp_A(A,B) -= slater_value_at_nucl(B,A)
       ! Projector
       do mu=1,mo_tot_num
-        cusp_A(A,B) += MOSlaOverlap_matrix(mu,B) * mo_value_at_nucl(mu,A)
+        cusp_A(A,B) += AO_orthoSlaOverlap_matrix(mu,B) * ao_ortho_value_at_nucl(mu,A)
       enddo
     enddo
    enddo
@@ -36,18 +36,28 @@ END_PROVIDER
 
    
 BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ]
-   implicit none
-   BEGIN_DOC
-   ! Equations to solve : A.C = B
-   END_DOC
-   
-   double precision, allocatable           :: AF(:,:)
-   integer :: info
-   allocate ( AF(nucl_num,nucl_num) )
-   
-   call get_pseudo_inverse(cusp_A,nucl_num,nucl_num,AF,size(AF,1))
-   call dgemm('N','N',nucl_num,mo_tot_num,nucl_num,1.d0, &
-     AF,size(AF,1), cusp_B, size(cusp_B,1), 0.d0, cusp_C, size(cusp_C,1))
+  implicit none
+  BEGIN_DOC
+  ! Equations to solve : A.C = B
+  END_DOC
+  
+  integer                        :: info
+  integer                        :: ipiv(nucl_num)
+  double precision, allocatable  :: AF(:,:)
+  allocate ( AF(nucl_num,nucl_num) )
+  
+  cusp_C(1:nucl_num,1:mo_tot_num) = cusp_B(1:nucl_num,1:mo_tot_num)
+  AF(1:nucl_num,1:nucl_num) = cusp_A(1:nucl_num,1:nucl_num)
+  call dgetrf(nucl_num,nucl_num,AF,size(AF,1),ipiv,info)
+  if (info /= 0) then
+    print *,  info
+    stop 'dgetrf failed'
+  endif
+  call dgetrs('N',nucl_num,mo_tot_num,AF,size(AF,1),ipiv,cusp_C,size(cusp_C,1),info)
+  if (info /= 0) then
+    print *,  info
+    stop 'dgetrs failed'
+  endif
 
 END_PROVIDER
 
diff --git a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f
index a8c1351a..8befdb91 100644
--- a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f
+++ b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f
@@ -76,21 +76,28 @@ subroutine run
   double precision               :: EHF
   integer                        :: i_it, i, j, k
    
-  mo_label = "CuspDressed"
+  mo_label = 'None'
 
 
-  print *,  HF_energy
-
-  call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(cusp_C)
+!  print *,  HF_energy
 
   do i=1,ao_num
     print *,  mo_coef(i,1), cusp_corrected_mos(i,1)
   enddo
   mo_coef(1:ao_num,1:mo_tot_num) = cusp_corrected_mos(1:ao_num,1:mo_tot_num)
-  TOUCH mo_coef
+  SOFT_TOUCH mo_coef slater_coef
+  call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(slater_coef)
+  call ezfio_set_Hartree_Fock_SlaterDressed_projector(ao_ortho_canonical_coef(1:ao_num,1:ao_num))
+  call ezfio_set_Hartree_Fock_SlaterDressed_ao_orthoSlaOverlap(AO_orthoSlaOverlap_matrix)
+  call save_mos
+  print *,  'ci'
+  print *, mo_coef(1:ao_num,1)
+  print *,  'cAi'
+  print *, slater_coef
 
-  EHF = HF_energy 
-  print *,  HF_energy
+
+!  EHF = HF_energy 
+!  print *,  HF_energy
 !  call Roothaan_Hall_SCF
   
 end
diff --git a/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f b/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f
index 910c5a16..afb85b64 100644
--- a/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f
+++ b/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f
@@ -28,14 +28,26 @@ BEGIN_PROVIDER [ double precision , ao_value_at_nucl, (ao_num,nucl_num) ]
   enddo
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, ao_ortho_value_at_nucl, (ao_num,nucl_num) ]
+  implicit none
+  BEGIN_DOC
+! Values of the molecular orbitals at the nucleus 
+  END_DOC
+
+  call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0,                    &
+      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),      &
+      ao_value_at_nucl, size(ao_value_at_nucl,1),                    &
+      0.d0, ao_ortho_value_at_nucl,size(ao_ortho_value_at_nucl,1))
+END_PROVIDER
+
 BEGIN_PROVIDER [ double precision, mo_value_at_nucl, (mo_tot_num,nucl_num) ]
   implicit none
   BEGIN_DOC
 ! Values of the molecular orbitals at the nucleus 
   END_DOC
 
-  call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0,                &
-      mo_coef_transp, size(mo_coef_transp,1),                        &
+  call dgemm('T','N',mo_tot_num,nucl_num,ao_num,1.d0,                &
+      mo_coef, size(mo_coef,1),                        &
       ao_value_at_nucl, size(ao_value_at_nucl,1),                    &
       0.d0, mo_value_at_nucl, size(mo_value_at_nucl,1))
 END_PROVIDER
@@ -54,10 +66,6 @@ BEGIN_PROVIDER [ double precision , slater_value_at_nucl, (nucl_num,nucl_num) ]
       x = nucl_coord(i,1) - nucl_coord(k,1)
       y = nucl_coord(i,2) - nucl_coord(k,2)
       z = nucl_coord(i,3) - nucl_coord(k,3)
-
-!      expo = slater_expo(i)*slater_expo(i)*((x*x) + (y*y) + (z*z))
-!      if (expo > 160.d0) cycle
-!      expo = dsqrt(expo)
       expo = slater_expo(i) * dsqrt((x*x) + (y*y) + (z*z))
       slater_value_at_nucl(i,k) = dexp(-expo) * slater_normalization(i)
     enddo
diff --git a/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f b/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f
new file mode 100644
index 00000000..f42a0b0f
--- /dev/null
+++ b/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f
@@ -0,0 +1,172 @@
+BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral_dressing, (ao_num_align,ao_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Dressing of the core hamiltonian in the orthogonal AO basis set
+  END_DOC
+
+  integer                        :: i,j,k
+  integer                        :: mu, nu, lambda, A
+  double precision               :: tmp
+
+  ao_ortho_mono_elec_integral_dressing = 0.d0
+  i = idx_dressing
+  do mu=1,ao_num
+    if (dabs(mo_coef_in_ao_ortho_basis(mu,i)) > 1.d-5) then
+      do A=1,nucl_num
+        tmp = 0.d0
+        do nu=1,ao_num
+          tmp += AO_orthoSlaOverlap_matrix(nu,A) * ao_ortho_mono_elec_integral(mu,nu)
+        enddo
+        ao_ortho_mono_elec_integral_dressing(mu,mu) += cusp_C(A,i) * (AO_orthoSlaH_matrix(mu,A) - tmp)
+      enddo
+      ao_ortho_mono_elec_integral_dressing(mu,mu) *= 1.d0/mo_coef_in_ao_ortho_basis(mu,i)
+    endif
+  enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral, (ao_num_align, ao_num) ]
+ implicit none
+ BEGIN_DOC
+ ! h core in orthogonal AO basis
+ END_DOC
+  double precision, allocatable ::  T(:,:)
+  allocate(T(ao_num,ao_num))
+  call dgemm('T','N',ao_num,ao_num,ao_num,1.d0,                    &
+      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),      &
+      ao_mono_elec_integral, size(ao_mono_elec_integral,1),          &
+      0.d0, T, size(T,1))
+  call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+      T, size(T,1),                                                  &
+      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),      &
+      0.d0, ao_ortho_mono_elec_integral, size(ao_ortho_mono_elec_integral,1))
+  deallocate(T)
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, ao_mono_elec_integral_dressing, (ao_num,ao_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Dressing of the core hamiltonian in the AO basis set
+  END_DOC
+  call ao_ortho_cano_to_ao(ao_ortho_mono_elec_integral_dressing,size(ao_ortho_mono_elec_integral_dressing,1),&
+    ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1))
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mo_mono_elec_integral_dressing, (mo_tot_num_align,mo_tot_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Dressing of the core hamiltonian in the MO basis set
+  END_DOC
+
+  call ao_to_mo(ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1),&
+    mo_mono_elec_integral_dressing,size(mo_mono_elec_integral_dressing,1))
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, idx_dressing ]
+ implicit none
+ BEGIN_DOC
+ ! Index of the MO which is being dressed
+ END_DOC
+ idx_dressing = 1
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, cusp_corrected_mos, (ao_num_align,mo_tot_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Dressing core hamiltonian in the AO basis set
+  END_DOC
+  integer :: i,j
+  double precision, allocatable :: F(:,:), M(:,:)
+  allocate(F(mo_tot_num_align,mo_tot_num),M(ao_num,mo_tot_num))
+
+  logical :: oneshot
+
+!  oneshot = .True. 
+  oneshot = .False.
+
+  if (oneshot) then
+    cusp_corrected_mos(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num)
+    slater_coef(1:nucl_num,1:mo_tot_num) = cusp_C(1:nucl_num,1:mo_tot_num)
+    return
+
+  else
+
+
+    do idx_dressing=1,mo_tot_num
+
+      if (idx_dressing>1) then
+        TOUCH idx_dressing
+      endif
+
+      do j=1,mo_tot_num
+        do i=1,mo_tot_num
+          F(i,j) = Fock_matrix_mo(i,j)
+        enddo
+      enddo
+
+      do j=1,mo_tot_num
+        do i=1,ao_num
+          M(i,j) = mo_coef(i,j)
+        enddo
+      enddo
+
+      integer :: it
+      do it=1,128
+
+  !      print *,  'X', ao_ortho_canonical_coef(1:ao_num,1:ao_num)
+  !      print *,  'C', mo_coef(1:ao_num,1:mo_tot_num)
+  !      print *,  'Cp', mo_coef_in_ao_ortho_basis(1:ao_num,1:mo_tot_num)
+  !      print *,  'cAi', cusp_C(1:nucl_num,1:mo_tot_num)
+  !      print *,  'FmuA', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num)
+  !      print *,  'Fock:', Fock_matrix_ao(1:ao_num,1:ao_num)
+  !      print *,  'Diag Dressing:', ao_ortho_mono_elec_integral_dressing(1:ao_num,1:ao_num)
+  !      print *,  'Dressing:', ao_mono_elec_integral_dressing(1:ao_num,1:ao_num)
+  !      print *,  'Dressed Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ao_mono_elec_integral_dressing(1:ao_num,1:ao_num)
+  !      print *,  'AO_orthoSlaOverlap_matrix', AO_orthoSlaOverlap_matrix(1:ao_num,1:nucl_num)
+  !      print *,  'AO_orthoSlaH_matrix', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num)
+  !      print *,  'ao_ortho_mono_elec_integral', ao_ortho_mono_elec_integral(1:ao_num,1:ao_num)
+  !      print *,  'Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num) 
+        do j=1,mo_tot_num
+          do i=1,mo_tot_num
+            Fock_matrix_mo(i,j) += mo_mono_elec_integral_dressing(i,j)
+          enddo
+        enddo
+        do i=1,mo_tot_num
+          Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i)
+        enddo
+  !      print *,  'Dressed Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num) 
+        double precision :: conv
+        conv = 0.d0
+        do j=1,mo_tot_num
+          do i=1,mo_tot_num
+            if (i==j) cycle
+            conv = max(conv,Fock_matrix_mo(i,j))
+          enddo
+        enddo
+        TOUCH Fock_matrix_mo Fock_matrix_diag_mo
+
+        mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_fock_matrix_mo(1:ao_num,1:mo_tot_num) 
+        TOUCH mo_coef
+  !print *,  'C', mo_coef(1:ao_num,1:mo_tot_num)
+  !print *,  '-----'
+        print *, idx_dressing, it, real(mo_coef(1,idx_dressing)), real(conv)
+        if (conv < 1.d-5) exit
+  !stop
+
+      enddo
+      cusp_corrected_mos(1:ao_num,idx_dressing) = mo_coef(1:ao_num,idx_dressing)
+      slater_coef(1:nucl_num,idx_dressing) = cusp_C(1:nucl_num,idx_dressing)
+    enddo
+
+    idx_dressing = 1
+    mo_coef(1:ao_num,1:mo_tot_num) = M(1:ao_num,1:mo_tot_num)
+    soft_TOUCH mo_coef idx_dressing slater_coef
+
+  endif
+
+END_PROVIDER
+
+
diff --git a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f
index 1111055b..680f3e1f 100644
--- a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f
+++ b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f
@@ -344,7 +344,7 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result)
   ss = k*ss
 
 ! Print result
-  write(*,*) ss
+!  write(*,*) ss
   result = 0.d0
 
 end
@@ -502,8 +502,8 @@ BEGIN_PROVIDER [ double precision, MOSla$X_matrix, (mo_tot_num, nucl_num) ]
   BEGIN_DOC
 ! <MO | Slater>
   END_DOC
-  call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0,                &
-      mo_coef_transp, size(mo_coef_transp,1),                        &
+  call dgemm('T','N',mo_tot_num,nucl_num,ao_num,1.d0,                &
+      mo_coef, size(mo_coef,1),                        &
       GauSla$X_matrix, size(GauSla$X_matrix,1),                    &
       0.d0, MOSla$X_matrix, size(MOSla$X_matrix,1))
 END_PROVIDER
@@ -520,6 +520,7 @@ BEGIN_PROVIDER [ double precision, AO_orthoSla$X_matrix, (ao_num, nucl_num) ]
  
 END_PROVIDER
 
+
 SUBST [ X ]
 
 Overlap ;;
@@ -616,9 +617,9 @@ BEGIN_PROVIDER [ double precision, AO_orthoSlaH_matrix, (ao_num, nucl_num) ]
   BEGIN_DOC
 ! <AO ortho | Slater>
   END_DOC
-  call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0,                &
-      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),                        &
-      GauSlaH_matrix, size(GauSlaH_matrix,1),                        &
-      0.d0, AO_orthoSlaH_matrix, size(AO_orthoSlaH_matrix,1))
+  call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0,                   &
+           ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), &
+           GauSlaH_matrix, size(GauSlaH_matrix,1),                   &
+           0.d0, AO_orthoSlaH_matrix, size(AO_orthoSlaH_matrix,1))
 END_PROVIDER
 
diff --git a/plugins/Hartree_Fock_SlaterDressed/slater.irp.f b/plugins/Hartree_Fock_SlaterDressed/slater.irp.f
index 8a11b9b1..27e4f10d 100644
--- a/plugins/Hartree_Fock_SlaterDressed/slater.irp.f
+++ b/plugins/Hartree_Fock_SlaterDressed/slater.irp.f
@@ -8,7 +8,10 @@ BEGIN_PROVIDER [ double precision, slater_expo, (nucl_num) ]
  if (exists) then
    slater_expo(1:nucl_num) = slater_expo_ezfio(1:nucl_num)
  else
-   slater_expo(1:nucl_num) = nucl_charge(1:nucl_num)
+   integer :: i
+   do i=1,nucl_num
+    slater_expo(i) = nucl_charge(i) 
+   enddo
    call ezfio_set_Hartree_Fock_SlaterDressed_slater_expo_ezfio(slater_expo)
  endif
 END_PROVIDER
diff --git a/src/Determinants/mo_energy_expval.irp.f b/src/Determinants/mo_energy_expval.irp.f
new file mode 100644
index 00000000..ea62a441
--- /dev/null
+++ b/src/Determinants/mo_energy_expval.irp.f
@@ -0,0 +1,159 @@
+BEGIN_PROVIDER [ double precision, mo_energy_expval, (N_states,mo_tot_num,2,2)]
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Third index is spin.
+  ! Fourth index is 1:creation, 2:annihilation
+  END_DOC
+  integer                        :: i,j,k
+  integer                        :: ispin, istate
+  integer                        :: hp
+  double precision               :: norm_out(N_states)
+  
+  integer, parameter             :: hole_particle(2) = (/ -1, 1 /)
+  double precision               :: energies(n_states)
+  
+  integer(bit_kind), allocatable :: psi_in_out(:,:,:)
+  double precision, allocatable  :: psi_in_out_coef(:,:)
+  double precision               :: E0(N_states), norm
+  double precision, parameter    :: t=1.d-3
+  
+  allocate (psi_in_out(N_int,2,N_det),psi_in_out_coef(N_det,N_states))
+  mo_energy_expval = 0.d0
+  
+  psi_in_out_coef(1:N_det,1:N_states) = psi_coef(1:N_det,1:N_states)
+  psi_in_out(1:N_int,1:2,1:N_det) = psi_det(1:N_int,1:2,1:N_det)
+  
+  ! Truncate the wave function
+  do istate=1,N_states
+    norm = 0.d0
+    do k=1,N_det
+      if (dabs(psi_in_out_coef(k,istate)) < t) then
+        psi_in_out_coef(k,istate) = 0.d0
+      endif
+      norm = norm + psi_in_out_coef(k,istate)*psi_in_out_coef(k,istate)
+    enddo
+    ASSERT (norm  > 0.d0)
+    norm = 1.d0/dsqrt(norm)
+    psi_in_out_coef(1:N_det,istate) = psi_in_out_coef(1:N_det,istate) * norm
+    call au0_h_au0(E0,psi_in_out,psi_in_out_coef,N_det,size(psi_in_out_coef,1))
+  enddo
+  
+  
+  do hp=1,2
+    do ispin=1,2
+      do i=1,mo_tot_num
+        psi_in_out_coef(1:N_det,1:N_states) = psi_coef(1:N_det,1:N_states)
+        psi_in_out(1:N_int,1:2,1:N_det) = psi_det(1:N_int,1:2,1:N_det)
+        call apply_exc_to_psi(i,hole_particle(hp),ispin,             &
+            norm_out,psi_in_out,psi_in_out_coef, N_det,N_det,N_det,N_states)
+        
+        ! Truncate the wave function
+        do istate=1,N_states
+          norm = 0.d0
+          do k=1,N_det
+            if (dabs(psi_in_out_coef(k,istate)) < t) then
+              psi_in_out_coef(k,istate) = 0.d0
+            endif
+            norm = norm + psi_in_out_coef(k,istate)*psi_in_out_coef(k,istate)
+          enddo
+          if (norm == 0.d0) then
+            cycle
+          endif
+          norm = 1.d0/dsqrt(norm)
+          psi_in_out_coef(1:N_det,istate) = psi_in_out_coef(1:N_det,istate) * norm
+        enddo
+        call au0_h_au0(energies,psi_in_out,psi_in_out_coef,N_det,size(psi_in_out_coef,1))
+        mo_energy_expval(1:N_states,i,ispin,hp) = energies(1:N_states) - E0(1:N_states)
+        print *,  i, ispin, real(energies(1)), real(E0(1))
+      enddo
+    enddo
+    
+  enddo
+  mo_energy_expval(1:N_states,1:mo_tot_num,1:2,1) = -mo_energy_expval(1:N_states,1:mo_tot_num,1:2,1) 
+
+END_PROVIDER
+
+  
+subroutine au0_h_au0(energies,psi_in,psi_in_coef,ndet,dim_psi_coef)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: ndet,dim_psi_coef
+  integer(bit_kind), intent(in)  :: psi_in(N_int,2,Ndet)
+  double precision,  intent(in)  :: psi_in_coef(dim_psi_coef,N_states)
+  double precision,  intent(out) :: energies(N_states)
+  
+  integer                        :: i,j, istate
+  double precision               :: hij,accu
+  double precision, allocatable  :: psi_coef_tmp(:)
+  
+  energies(1:N_states) = 0.d0
+  do i = 1, Ndet
+    if(sum(dabs(psi_in_coef(i,1:N_states)))==0.d0) then
+      cycle
+    endif
+    call diag_H_mat_elem_au0_h_au0(psi_in(1,1,i),N_int,hij)
+    do istate=1,N_states
+      energies(istate) += psi_in_coef(i,istate) * psi_in_coef(i,istate) * hij
+    enddo
+    do j = i+1, Ndet
+      if(sum(dabs(psi_in_coef(j,1:N_states)))==0.d0) then
+        cycle
+      endif
+      call i_H_j(psi_in(1,1,i),psi_in(1,1,j),N_int,hij)
+      hij = hij+hij
+      do istate=1,N_states
+        energies(istate) = energies(istate) + psi_in_coef(i,istate) * psi_in_coef(j,istate) * hij
+      enddo
+    enddo
+  enddo
+end
+
+subroutine diag_H_mat_elem_au0_h_au0(det_in,Nint,hii)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Computes <i|H|i> for any determinant i. Used for wave functions with an additional electron.
+  END_DOC
+  integer,intent(in)             :: Nint
+  integer(bit_kind),intent(in)   :: det_in(Nint,2)
+  double precision, intent(out)  :: hii
+  
+  integer                        :: i, j, iorb, jorb
+  integer                        :: occ(Nint*bit_kind_size,2)
+  integer                        :: elec_num_tab_local(2)
+  
+  hii = 0.d0
+  call bitstring_to_list(det_in(1,1), occ(1,1), elec_num_tab_local(1), Nint)
+  call bitstring_to_list(det_in(1,2), occ(1,2), elec_num_tab_local(2), Nint)
+  
+  ! alpha - alpha
+  do i = 1, elec_num_tab_local(1)
+    iorb =  occ(i,1)
+    hii += mo_mono_elec_integral(iorb,iorb)
+    do j = i+1, elec_num_tab_local(1)
+      jorb = occ(j,1)
+      hii +=  mo_bielec_integral_jj_anti(jorb,iorb)
+    enddo
+  enddo
+  
+  ! beta - beta
+  do i = 1, elec_num_tab_local(2)
+    iorb =  occ(i,2)
+    hii += mo_mono_elec_integral(iorb,iorb)
+    do j = i+1, elec_num_tab_local(2)
+      jorb = occ(j,2)
+      hii +=  mo_bielec_integral_jj_anti(jorb,iorb)
+    enddo
+  enddo
+  
+  ! alpha - beta
+  do i = 1, elec_num_tab_local(2)
+    iorb =  occ(i,2)
+    do j = 1, elec_num_tab_local(1)
+      jorb = occ(j,1)
+      hii +=  mo_bielec_integral_jj(jorb,iorb)
+    enddo
+  enddo
+  
+end
diff --git a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
index 644f89c2..f24e4974 100644
--- a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
+++ b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
@@ -56,7 +56,6 @@ subroutine push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value,
     stop 'error'
   endif
 
-! Activate is zmq_socket_push is a REQ
 IRP_IF ZMQ_PUSH
 IRP_ELSE
   integer :: idummy
@@ -189,7 +188,6 @@ subroutine ao_bielec_integrals_in_map_collector
 
       rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
 
-! Activate if zmq_socket_pull is a REP
 IRP_IF ZMQ_PUSH
 IRP_ELSE
       rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)
diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f
index b0eabfbd..d77374a8 100644
--- a/src/MO_Basis/ao_ortho_canonical.irp.f
+++ b/src/MO_Basis/ao_ortho_canonical.irp.f
@@ -76,19 +76,20 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_
  ! AO_cart_to_sphe_coef^(-1)
  END_DOC
 
- call get_pseudo_inverse(ao_cart_to_sphe_coef,ao_num,ao_cart_to_sphe_num, &
-   ao_cart_to_sphe_inv, size(ao_cart_to_sphe_coef,1))
+ call get_pseudo_inverse(ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1),&
+   ao_num,ao_cart_to_sphe_num, &
+   ao_cart_to_sphe_inv, size(ao_cart_to_sphe_inv,1))
 END_PROVIDER
 
 
 
-BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef_inv, (ao_num,ao_num)]
+BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef_inv, (ao_num_align,ao_num)]
  implicit none
  BEGIN_DOC
 ! ao_ortho_canonical_coef^(-1)
  END_DOC
- call get_pseudo_inverse(ao_ortho_canonical_coef,ao_num,ao_num, &
-   ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef,1))
+ call get_inverse(ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),&
+     ao_num, ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1))
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)]
@@ -100,11 +101,15 @@ END_PROVIDER
 ! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
   END_DOC
   integer :: i
-  ao_ortho_canonical_coef(:,:) = 0.d0
+  ao_ortho_canonical_coef = 0.d0
   do i=1,ao_num
     ao_ortho_canonical_coef(i,i) = 1.d0
   enddo
 
+call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num)
+ao_ortho_canonical_num=ao_num
+return
+
   if (ao_cartesian) then
 
     ao_ortho_canonical_num = ao_num
diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f
index 260af392..8962dd00 100644
--- a/src/MO_Basis/mos.irp.f
+++ b/src/MO_Basis/mos.irp.f
@@ -72,8 +72,10 @@ BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num_align, mo_
  implicit none
  BEGIN_DOC
  ! MO coefficients in orthogonalized AO basis
+ !
+ ! C^(-1).C_mo
  END_DOC
- call dgemm('T','N',ao_num,mo_tot_num,ao_num,1.d0,                   &
+ call dgemm('N','N',ao_num,mo_tot_num,ao_num,1.d0,                   &
      ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),&
      mo_coef, size(mo_coef,1), 0.d0,                                 &
      mo_coef_in_ao_ortho_basis, size(mo_coef_in_ao_ortho_basis,1))
@@ -155,6 +157,8 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
   implicit none
   BEGIN_DOC
   ! Transform A from the AO basis to the MO basis
+  !
+  ! C.A_ao.Ct
   END_DOC
   integer, intent(in)            :: LDA_ao,LDA_mo
   double precision, intent(in)   :: A_ao(LDA_ao)
@@ -181,6 +185,8 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
   implicit none
   BEGIN_DOC
   ! Transform A from the MO basis to the AO basis
+  !
+  ! (S.C).A_mo.(S.C)t
   END_DOC
   integer, intent(in)            :: LDA_ao,LDA_mo
   double precision, intent(in)   :: A_mo(LDA_mo)
@@ -273,6 +279,8 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)
   implicit none
   BEGIN_DOC
   ! Transform A from the AO basis to the orthogonal AO basis
+  !
+  ! C^(-1).A_ao.Ct^(-1)
   END_DOC
   integer, intent(in)            :: LDA_ao,LDA
   double precision, intent(in)   :: A_ao(LDA_ao,*)
@@ -296,35 +304,3 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)
   deallocate(T)
 end
 
-subroutine mo_to_ao_ortho_cano(A_mo,LDA_mo,A_ao,LDA_ao)
-  implicit none
-  BEGIN_DOC
-  ! Transform A from the AO orthogonal basis to the AO basis
-  END_DOC
-  integer, intent(in)            :: LDA_ao,LDA_mo
-  double precision, intent(in)   :: A_mo(LDA_mo)
-  double precision, intent(out)  :: A_ao(LDA_ao)
-  double precision, allocatable  :: T(:,:), SC(:,:)
-  
-  allocate ( SC(ao_num_align,mo_tot_num) )
-  allocate ( T(mo_tot_num_align,ao_num) )
-  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
-  
-  call dgemm('N','N', ao_num, mo_tot_num, ao_num,                    &
-      1.d0, ao_overlap,size(ao_overlap,1),                           &
-      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),                                      &
-      0.d0, SC, ao_num_align)
-  
-  call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num,                &
-      1.d0, A_mo,LDA_mo,                                             &
-      SC, size(SC,1),                                      &
-      0.d0, T, mo_tot_num_align)
-  
-  call dgemm('N','N', ao_num, ao_num, mo_tot_num,                    &
-      1.d0, SC,size(SC,1),                                 &
-      T, mo_tot_num_align,                                           &
-      0.d0, A_ao, LDA_ao)
-  
-  deallocate(T,SC)
-end
-
diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f
index 92bcdf53..f09e98ff 100644
--- a/src/Utils/LinearAlgebra.irp.f
+++ b/src/Utils/LinearAlgebra.irp.f
@@ -72,8 +72,6 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
   double precision, allocatable  :: S_half(:,:)
   !DEC$ ATTRIBUTES ALIGN : 64    :: U, Vt, D
   integer                        :: info, i, j
-!call ortho_lowdin(overlap,LDA,N,C,LDC,m)
-!return
   
   if (n < 2) then
     return
@@ -270,14 +268,42 @@ end
 
 
 
-subroutine get_pseudo_inverse(A,m,n,C,LDA)
+subroutine get_inverse(A,LDA,m,C,LDC)
+  implicit none
+  BEGIN_DOC
+  ! Returns the inverse of the square matrix A
+  END_DOC
+  integer, intent(in)            :: m, LDA, LDC
+  double precision, intent(in)   :: A(LDA,m)
+  double precision, intent(out)  :: C(LDC,m)
+
+  integer                        :: info,lwork
+  integer, allocatable           :: ipiv(:)
+  double precision,allocatable   :: work(:)
+  allocate (ipiv(ao_num), work(ao_num*ao_num))
+  lwork = size(work)
+  C(1:m,1:m) = A(1:m,1:m)
+  call dgetrf(m,m,C,size(C,1),ipiv,info)
+  if (info /= 0) then
+    print *,  info
+    stop 'error in inverse (dgetrf)'
+  endif
+  call dgetri(m,C,size(C,1),ipiv,work,lwork,info)
+  if (info /= 0) then
+    print *,  info
+    stop 'error in inverse (dgetri)'
+  endif
+  deallocate(ipiv,work)
+end
+
+subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC)
   implicit none
   BEGIN_DOC
   ! Find C = A^-1
   END_DOC
-  integer, intent(in)            :: m,n, LDA
+  integer, intent(in)            :: m,n, LDA, LDC
   double precision, intent(in)   :: A(LDA,n)
-  double precision, intent(out)  :: C(n,m)
+  double precision, intent(out)  :: C(LDC,m)
   
   double precision, allocatable  :: U(:,:), D(:), Vt(:,:), work(:), A_tmp(:,:)
   integer                        :: info, lwork
@@ -304,7 +330,7 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA)
   endif
   
   do i=1,n
-    if (dabs(D(i)) > 1.d-6) then
+    if (D(i)/D(1) > 1.d-10) then
       D(i) = 1.d0/D(i)
     else
       D(i) = 0.d0
@@ -315,7 +341,7 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA)
   do i=1,m
     do j=1,n
       do k=1,n
-        C(j,i) += U(i,k) * D(k) * Vt(k,j)
+        C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j)
       enddo
     enddo
   enddo
diff --git a/src/ZMQ/utils.irp.f b/src/ZMQ/utils.irp.f
index 4b0cd0c5..7b26ed57 100644
--- a/src/ZMQ/utils.irp.f
+++ b/src/ZMQ/utils.irp.f
@@ -233,7 +233,7 @@ function new_zmq_pull_socket()
      stop 'zmq_context is uninitialized'
   endif
 IRP_IF ZMQ_PUSH
- new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL)
+  new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL)
 IRP_ELSE
   new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_REP)
 IRP_ENDIF