diff --git a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f
index 35d9bcc4..736d1a97 100644
--- a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f
+++ b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f
@@ -1,7 +1,7 @@
 BEGIN_PROVIDER [double precision, HF_density_matrix_ao_alpha, (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! Alpha density matrix in the AO basis 
+   ! S^{-1}.P_alpha.S^{-1}
    END_DOC
    
    call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, &
@@ -14,7 +14,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! Beta density matrix in the AO basis
+   ! S^{-1}.P_beta.S^{-1}
    END_DOC
    
    call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, &
@@ -27,7 +27,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! Density matrix in the AO basis
+   ! S^{-1}.P.S^{-1}  where P = C.C^t
    END_DOC
    ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1))
    if (elec_alpha_num== elec_beta_num) then
diff --git a/plugins/Hartree_Fock/huckel.irp.f b/plugins/Hartree_Fock/huckel.irp.f
index 1ecf7412..c9e32ad5 100644
--- a/plugins/Hartree_Fock/huckel.irp.f
+++ b/plugins/Hartree_Fock/huckel.irp.f
@@ -7,11 +7,11 @@ subroutine huckel_guess
   double precision               :: accu
   double precision               :: c
   character*(64)                 :: label
-  double precision, allocatable  :: A(:,:), T(:,:)
+  double precision, allocatable  :: A(:,:)
   label = "Guess"
-  c = 0.25d0 * 1.75d0
+  c = 0.5d0 * 1.75d0
 
-  allocate (A(ao_num, ao_num),T(ao_num, ao_num))
+  allocate (A(ao_num, ao_num))
   A = 0.d0
   do j=1,ao_num
     do i=1,ao_num
@@ -29,7 +29,6 @@ subroutine huckel_guess
   mo_coef = eigenvectors_fock_matrix_mo
   SOFT_TOUCH mo_coef
   call save_mos
-  print *,  'E=', HF_energy
-  deallocate(A,T)
+  deallocate(A)
 
 end
diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f
index fd3fa8ec..ecb149b9 100644
--- a/src/AO_Basis/ao_overlap.irp.f
+++ b/src/AO_Basis/ao_overlap.irp.f
@@ -34,6 +34,8 @@
     power_A(1)  = ao_power( j, 1 )
     power_A(2)  = ao_power( j, 2 )
     power_A(3)  = ao_power( j, 3 )
+    !DEC$ VECTOR ALIGNED
+    !DEC$ VECTOR ALWAYS
     do i= 1,ao_num
       ao_overlap(i,j)= 0.d0
       ao_overlap_x(i,j)= 0.d0
@@ -47,6 +49,7 @@
       power_B(3)  = ao_power( i, 3 )
       do n = 1,ao_prim_num(j)
       alpha = ao_expo_ordered_transp(n,j)
+      !DEC$ VECTOR ALIGNED
       do l = 1, ao_prim_num(i)
         beta = ao_expo_ordered_transp(l,i)
         call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
@@ -100,6 +103,8 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
    power_A(1)  = ao_power( j, 1 )
    power_A(2)  = ao_power( j, 2 )
    power_A(3)  = ao_power( j, 3 )
+   !DEC$ VECTOR ALIGNED
+   !DEC$ VECTOR ALWAYS
    do i= 1,ao_num
     ao_overlap_abs(i,j)= 0.d0
     B_center(1) = nucl_coord( ao_nucl(i), 1 )
@@ -110,6 +115,7 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
     power_B(3)  = ao_power( i, 3 )
     do n = 1,ao_prim_num(j)
      alpha = ao_expo_ordered_transp(n,j)
+     !DEC$ VECTOR ALIGNED
      do l = 1, ao_prim_num(i)
       beta = ao_expo_ordered_transp(l,i)
       call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1)
diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg
index 6bf6faff..f4c5d866 100644
--- a/src/Determinants/EZFIO.cfg
+++ b/src/Determinants/EZFIO.cfg
@@ -50,7 +50,7 @@ default: 0.99
 type: Threshold
 doc: Thresholds on selectors (fraction of the norm) 
 interface: ezfio,provider,ocaml
-default:        0.995
+default:        0.999
 
 [n_int]
 interface: ezfio
diff --git a/src/Integrals_Monoelec/ao_mono_ints.irp.f b/src/Integrals_Monoelec/ao_mono_ints.irp.f
index 8cfd25cf..4646326e 100644
--- a/src/Integrals_Monoelec/ao_mono_ints.irp.f
+++ b/src/Integrals_Monoelec/ao_mono_ints.irp.f
@@ -7,6 +7,7 @@
  ! : sum of the kinetic and nuclear electronic potential 
   END_DOC
   do j = 1, ao_num
+   !DIR$ VECTOR ALIGNED
    do i = 1, ao_num
     ao_mono_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + ao_kinetic_integral(i,j) + ao_pseudo_integral(i,j)
    enddo
diff --git a/src/Integrals_Monoelec/kin_ao_ints.irp.f b/src/Integrals_Monoelec/kin_ao_ints.irp.f
index da4de4d7..6cb2aa49 100644
--- a/src/Integrals_Monoelec/kin_ao_ints.irp.f
+++ b/src/Integrals_Monoelec/kin_ao_ints.irp.f
@@ -45,6 +45,8 @@
    power_A(1)  = ao_power( j, 1 )
    power_A(2)  = ao_power( j, 2 )
    power_A(3)  = ao_power( j, 3 )
+   !DEC$ VECTOR ALIGNED
+   !DEC$ VECTOR ALWAYS
    do i= 1,ao_num
     ao_deriv2_x(i,j)= 0.d0
     ao_deriv2_y(i,j)= 0.d0
@@ -57,6 +59,7 @@
     power_B(3)  = ao_power( i, 3 )
     do n = 1,ao_prim_num(j)
      alpha = ao_expo_ordered_transp(n,j)
+     !DEC$ VECTOR ALIGNED
      do l = 1, ao_prim_num(i)
       beta = ao_expo_ordered_transp(l,i)
       call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1)
@@ -136,6 +139,8 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)]
     !$OMP  PRIVATE(i,j) &
     !$OMP  SHARED(ao_num, ao_num_align, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z)
     do j = 1, ao_num
+      !DEC$ VECTOR ALWAYS
+      !DEC$ VECTOR ALIGNED
       do i = 1, ao_num
       ao_kinetic_integral(i,j) = -0.5d0 * (ao_deriv2_x(i,j) + ao_deriv2_y(i,j) + ao_deriv2_z(i,j) )
       enddo
diff --git a/src/Integrals_Monoelec/spread_dipole_ao.irp.f b/src/Integrals_Monoelec/spread_dipole_ao.irp.f
index c510d58b..5611ec7f 100644
--- a/src/Integrals_Monoelec/spread_dipole_ao.irp.f
+++ b/src/Integrals_Monoelec/spread_dipole_ao.irp.f
@@ -35,6 +35,8 @@
    power_A(1)  = ao_power( j, 1 )
    power_A(2)  = ao_power( j, 2 )
    power_A(3)  = ao_power( j, 3 )
+   !DEC$ VECTOR ALIGNED
+   !DEC$ VECTOR ALWAYS
    do i= 1,ao_num
     B_center(1) = nucl_coord( ao_nucl(i), 1 )
     B_center(2) = nucl_coord( ao_nucl(i), 2 )
@@ -47,6 +49,7 @@
     accu_z = 0.d0
     do n = 1,ao_prim_num(j)
      alpha = ao_expo_ordered_transp(n,j)
+     !DEC$ VECTOR ALIGNED
      do l = 1, ao_prim_num(i)
       c = ao_coef_normalized_ordered_transp(n,j)*ao_coef_normalized_ordered_transp(l,i)
       beta = ao_expo_ordered_transp(l,i)
@@ -106,6 +109,8 @@
    power_A(1)  = ao_power( j, 1 )
    power_A(2)  = ao_power( j, 2 )
    power_A(3)  = ao_power( j, 3 )
+   !DEC$ VECTOR ALIGNED
+   !DEC$ VECTOR ALWAYS
    do i= 1,ao_num
     B_center(1) = nucl_coord( ao_nucl(i), 1 )
     B_center(2) = nucl_coord( ao_nucl(i), 2 )
@@ -118,6 +123,7 @@
     accu_z = 0.d0
     do n = 1,ao_prim_num(j)
      alpha = ao_expo_ordered_transp(n,j)
+     !DEC$ VECTOR ALIGNED
      do l = 1, ao_prim_num(i)
       beta = ao_expo_ordered_transp(l,i)
       c = ao_coef_normalized_ordered_transp(l,i)*ao_coef_normalized_ordered_transp(n,j)
@@ -177,6 +183,8 @@
    power_A(1)  = ao_power( j, 1 )
    power_A(2)  = ao_power( j, 2 )
    power_A(3)  = ao_power( j, 3 )
+   !DEC$ VECTOR ALIGNED
+   !DEC$ VECTOR ALWAYS
    do i= 1,ao_num
     B_center(1) = nucl_coord( ao_nucl(i), 1 )
     B_center(2) = nucl_coord( ao_nucl(i), 2 )
@@ -189,6 +197,7 @@
     accu_z = 0.d0
     do n = 1,ao_prim_num(j)
      alpha = ao_expo_ordered_transp(n,j)
+     !DEC$ VECTOR ALIGNED
      do l = 1, ao_prim_num(i)
       beta = ao_expo_ordered_transp(l,i)
       call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f
index c2cbc5f2..4bcd7221 100644
--- a/src/MO_Basis/mos.irp.f
+++ b/src/MO_Basis/mos.irp.f
@@ -49,6 +49,9 @@ END_PROVIDER
       do j=1,ao_num
         mo_coef(j,i) = buffer(j,i)
       enddo
+      do j=ao_num+1,ao_num_align
+        mo_coef(j,i) = 0.d0
+      enddo
     enddo
     deallocate(buffer)
   else
@@ -57,6 +60,9 @@ END_PROVIDER
       do j=1,ao_num
         mo_coef(j,i) = ao_ortho_canonical_coef(j,i)
       enddo
+      do j=ao_num+1,ao_num_align
+        mo_coef(j,i) = 0.d0
+      enddo
     enddo
   endif
 END_PROVIDER
@@ -202,6 +208,32 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
   deallocate(T)
 end
 
+subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
+  implicit none
+  BEGIN_DOC
+  ! Transform A from the MO basis to the S^-1 AO basis
+  ! Useful for density matrix
+  END_DOC
+  integer, intent(in)            :: LDA_ao,LDA_mo
+  double precision, intent(in)   :: A_mo(LDA_mo,mo_tot_num)
+  double precision, intent(out)  :: A_ao(LDA_ao,ao_num)
+  double precision, allocatable  :: T(:,:)
+  
+  allocate ( T(mo_tot_num_align,ao_num) )
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
+  
+  call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num,                &
+      1.d0, A_mo,size(A_mo,1),                                       &
+      mo_coef, size(mo_coef,1),                                      &
+      0.d0, T, size(T,1))
+  
+  call dgemm('N','N', ao_num, ao_num, mo_tot_num,                    &
+      1.d0, mo_coef,size(mo_coef,1),                                 &
+      T, size(T,1),                                                  &
+      0.d0, A_ao, size(A_ao,1))
+  
+  deallocate(T)
+end
 
 subroutine mix_mo_jk(j,k)
  implicit none