diff --git a/plugins/Hartree_Fock/diagonalize_fock.irp.f b/plugins/Hartree_Fock/diagonalize_fock.irp.f
index 2aabf9bd..78d08aa3 100644
--- a/plugins/Hartree_Fock/diagonalize_fock.irp.f
+++ b/plugins/Hartree_Fock/diagonalize_fock.irp.f
@@ -1,38 +1,46 @@
- BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (ao_num) ]
+ BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (mo_tot_num) ]
 &BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num_align,mo_tot_num) ]
    implicit none
    BEGIN_DOC
    ! Diagonal Fock matrix in the MO basis
    END_DOC
    
-   integer                        :: i,j
+   integer                        :: i,j, m
    integer                        :: liwork, lwork, n, info
-   integer, allocatable           :: iwork(:)
-   double precision, allocatable  :: work(:), F(:,:), S(:,:)
+   integer, allocatable           :: iwork(:), isuppz(:)
+   double precision, allocatable  :: work(:), F(:,:), F2(:,:)
+   integer                        :: iorb,jorb
    
    
-   allocate( F(mo_tot_num,mo_tot_num) )
+   allocate( F(mo_tot_num,mo_tot_num),F2(mo_tot_num,mo_tot_num), isuppz(2*mo_tot_num) )
    do j=1,mo_tot_num
      do i=1,mo_tot_num
        F(i,j) = Fock_matrix_mo(i,j)
      enddo
    enddo
    if(no_oa_or_av_opt)then
-     integer                        :: iorb,jorb
      do i = 1, n_act_orb
        iorb = list_act(i)
+       ASSERT (iorb > 0)
+       ASSERT (iorb <= mo_tot_num)
        do j = 1, n_inact_orb
          jorb = list_inact(j)
+        ASSERT (jorb > 0)
+        ASSERT (jorb <= mo_tot_num)
          F(iorb,jorb) = 0.d0
          F(jorb,iorb) = 0.d0
        enddo
        do j = 1, n_virt_orb
          jorb = list_virt(j)
+          ASSERT (jorb > 0)
+          ASSERT (jorb <= mo_tot_num)
          F(iorb,jorb) = 0.d0
          F(jorb,iorb) = 0.d0
        enddo
        do j = 1, n_core_orb
          jorb = list_core(j)
+          ASSERT (jorb > 0)
+          ASSERT (jorb <= mo_tot_num)
          F(iorb,jorb) = 0.d0
          F(jorb,iorb) = 0.d0
        enddo
@@ -53,57 +61,28 @@
   
   n = mo_tot_num
   lwork = 1+6*n + 2*n*n
-  liwork = 3 + 5*n
+  liwork = 10*n
   
   allocate(work(lwork))
   allocate(iwork(liwork) )
   
-  lwork = -1
-  liwork = -1
-  
-  
-  integer                        :: m, ISUPPZ(mo_tot_num)
-  call dsyevr('V', 'A', 'U', mo_tot_num, F, size(F,1),                &
-      -100.d0, 100.d0, 1, mo_tot_num, 0.d0,                           &
-      m, diagonal_Fock_matrix_mo,                                     &
-      eigenvectors_Fock_matrix_mo,                                    &
-      size(eigenvectors_Fock_matrix_mo,1),                            &
+  call dsyevr('V', 'A', 'U', mo_tot_num, F, size(F,1),               &
+      -100.d0, 100.d0, 1, mo_tot_num, 0.d0,                          &
+      m, diagonal_Fock_matrix_mo,                                    &
+      F2, size(F2,1),                                                &
       isuppz, work, lwork, iwork, liwork, info)
   
-  
-  
-  if (info /= 0) then
-    print *,  irp_here//' DSYEVD failed : ', info
-    stop 1
-  endif
-  lwork = int(work(1))
-  liwork = iwork(1)
-  deallocate(iwork)
-  deallocate(work)
-  
-  allocate(work(lwork))
-  allocate(iwork(liwork) )
-  
-  call dsyevr('V', 'A', 'U', mo_tot_num, F, size(F,1),                &
-      -100.d0, 100.d0, 1, mo_tot_num, 0.d0,                           &
-      m, diagonal_Fock_matrix_mo,                                     &
-      eigenvectors_Fock_matrix_mo,                                    &
-      size(eigenvectors_Fock_matrix_mo,1),                            &
-      isuppz, work, lwork, iwork, liwork, info)
-  
-  deallocate(iwork)
-  
   if (info /= 0) then
     print *,  irp_here//' DSYEV failed : ', info
     stop 1
   endif
   
-  F(1:mo_tot_num,1:mo_tot_num) = eigenvectors_Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num)
-  
-  call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0,              &
-      mo_coef, size(mo_coef,1), F, size(F,1),                         &
+  call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0,             &
+      mo_coef, size(mo_coef,1), F2, size(F2,1),                      &
       0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1))
-  deallocate(work, F)
+  deallocate(work, F2, F)
+  deallocate(iwork, isuppz)
+