diff --git a/src/BiInts/ao_bi_integrals.irp.f b/src/BiInts/ao_bi_integrals.irp.f
index 904122d6..fbd6fb98 100644
--- a/src/BiInts/ao_bi_integrals.irp.f
+++ b/src/BiInts/ao_bi_integrals.irp.f
@@ -194,6 +194,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
   use map_module
   BEGIN_DOC
   !  Map of Atomic integrals
+  !     i(r1) j(r2) 1/r12 k(r1) l(r2)
   END_DOC
   
   integer                        :: i,j,k,l
diff --git a/src/Hartree_Fock/Fock_matrix_mo.irp.f b/src/Hartree_Fock/Fock_matrix_mo.irp.f
index 38b19362..d9e05510 100644
--- a/src/Hartree_Fock/Fock_matrix_mo.irp.f
+++ b/src/Hartree_Fock/Fock_matrix_mo.irp.f
@@ -1,4 +1,4 @@
- BEGIN_PROVIDER [ double precision, Fock_matrix_mo, (mo_tot_num,mo_tot_num) ]
+ BEGIN_PROVIDER [ double precision, Fock_matrix_mo, (mo_tot_num_align,mo_tot_num) ]
 &BEGIN_PROVIDER [ double precision, Fock_matrix_diag_mo, (mo_tot_num)]
    implicit none
    BEGIN_DOC
@@ -17,7 +17,6 @@
    !
    END_DOC
    integer                        :: i,j,n
-   double precision               :: get_mo_bielec_integral
    if (elec_alpha_num == elec_beta_num) then
      Fock_matrix_mo = Fock_matrix_alpha_mo
    else
diff --git a/src/Hartree_Fock/HF_density_matrix_ao.irp.f b/src/Hartree_Fock/HF_density_matrix_ao.irp.f
new file mode 100644
index 00000000..b0d58344
--- /dev/null
+++ b/src/Hartree_Fock/HF_density_matrix_ao.irp.f
@@ -0,0 +1,87 @@
+ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ]
+&BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num_align,ao_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Alpha and Beta density matrix in the AO basis
+   END_DOC
+   integer                        :: i,j,k,l1,l2
+   integer, allocatable           :: mo_occ(:,:)
+   
+   allocate ( mo_occ(elec_alpha_num,2) )
+   call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int)
+   ASSERT ( j==elec_alpha_num )
+   
+   call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int)
+   ASSERT ( j==elec_beta_num )
+   
+   do j=1,ao_num
+     !DIR$ VECTOR ALIGNED
+     do i=1,ao_num_align
+       HF_density_matrix_ao_alpha(i,j) = 0.d0
+       HF_density_matrix_ao_beta (i,j) = 0.d0
+     enddo
+     do k=1,elec_beta_num
+       l1 = mo_occ(k,1)
+       l2 = mo_occ(k,2)
+       !DIR$ VECTOR ALIGNED
+       do i=1,ao_num
+         HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
+             mo_coef(i,l1) * mo_coef(j,l1)
+         HF_density_matrix_ao_beta (i,j) = HF_density_matrix_ao_beta (i,j) +&
+             mo_coef(i,l2) * mo_coef(j,l2)
+       enddo
+     enddo
+     do k=elec_beta_num+1,elec_alpha_num
+       l1 = mo_occ(k,1)
+       !DIR$ VECTOR ALIGNED
+       do i=1,ao_num
+         HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
+             mo_coef(i,l1) * mo_coef(j,l1)
+       enddo
+     enddo
+   enddo
+   deallocate(mo_occ)
+END_PROVIDER
+ 
+BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Density matrix in the AO basis
+   END_DOC
+   integer                        :: i,j,k,l1,l2
+   integer, allocatable           :: mo_occ(:,:)
+   
+   allocate ( mo_occ(elec_alpha_num,2) )
+   call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int)
+   ASSERT ( j==elec_alpha_num )
+   
+   call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int)
+   ASSERT ( j==elec_beta_num )
+   
+   do j=1,ao_num
+     !DIR$ VECTOR ALIGNED
+     do i=1,ao_num_align
+       HF_density_matrix_ao(i,j) = 0.d0
+     enddo
+     do k=1,elec_beta_num
+       l1 = mo_occ(k,1)
+       l2 = mo_occ(k,2)
+       !DIR$ VECTOR ALIGNED
+       do i=1,ao_num
+         HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) +           &
+             mo_coef(i,l1) * mo_coef(j,l1) +                         &
+             mo_coef(i,l2) * mo_coef(j,l2)
+       enddo
+     enddo
+     do k=elec_beta_num+1,elec_alpha_num
+       l1 = mo_occ(k,1)
+       !DIR$ VECTOR ALIGNED
+       do i=1,ao_num
+         HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) +           &
+             mo_coef(i,l1) * mo_coef(j,l1)
+       enddo
+     enddo
+   enddo
+   deallocate(mo_occ)
+END_PROVIDER
+ 
diff --git a/src/Hartree_Fock/hartree_fock.ezfio_config b/src/Hartree_Fock/hartree_fock.ezfio_config
index bbb383f9..134486f4 100644
--- a/src/Hartree_Fock/hartree_fock.ezfio_config
+++ b/src/Hartree_Fock/hartree_fock.ezfio_config
@@ -1,4 +1,6 @@
 hartree_fock
   thresh_scf       double precision
   n_it_scf_max     integer
+  direct           logical
+  diis             logical
 
diff --git a/src/Hartree_Fock/options.irp.f b/src/Hartree_Fock/options.irp.f
index b3770418..303ce89d 100644
--- a/src/Hartree_Fock/options.irp.f
+++ b/src/Hartree_Fock/options.irp.f
@@ -16,7 +16,7 @@ BEGIN_PROVIDER [ double precision,thresh_SCF ]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ integer ,n_it_scf_max]
+BEGIN_PROVIDER [ integer, n_it_scf_max]
   implicit none
   BEGIN_DOC  
 !  Maximum number of SCF iterations
@@ -34,3 +34,40 @@ BEGIN_PROVIDER [ integer ,n_it_scf_max]
 
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ logical, do_direct_SCF ]
+  implicit none
+  BEGIN_DOC  
+! If True, compute integrals on the fly
+  END_DOC
+
+  logical :: has
+  PROVIDE ezfio_filename
+  call ezfio_has_Hartree_Fock_direct(has)
+  if (has) then
+    call ezfio_get_Hartree_Fock_direct(do_direct_SCF)
+  else
+    do_direct_SCF = .False.
+    call ezfio_set_Hartree_Fock_direct(do_direct_SCF)
+  endif
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ logical, do_DIIS ]
+  implicit none
+  BEGIN_DOC  
+! If True, compute integrals on the fly
+  END_DOC
+
+  logical :: has
+  PROVIDE ezfio_filename
+  call ezfio_has_Hartree_Fock_DIIS(has)
+  if (has) then
+    call ezfio_get_Hartree_Fock_DIIS(do_DIIS)
+  else
+    do_DIIS = .False.
+    call ezfio_set_Hartree_Fock_DIIS(do_DIIS)
+  endif
+
+END_PROVIDER
+
diff --git a/src/Makefile.common b/src/Makefile.common
index 4c2b2b2b..2dad85e2 100644
--- a/src/Makefile.common
+++ b/src/Makefile.common
@@ -113,7 +113,7 @@ clean_links:
 endif
 
 LIB+=$(EZFIO) $(MKL)
-IRPF90+=$(patsubst %, -I %, $(INCLUDE_DIRS)) 
+IRPF90+=$(patsubst %, -I %, $(INCLUDE_DIRS)) $(IRPF90_FLAGS)
 
 irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile $(EZFIO) NEEDED_MODULES
 	$(IRPF90)
diff --git a/src/MonoInts/ao_mono_ints.irp.f b/src/MonoInts/ao_mono_ints.irp.f
index f77924d3..9a135dda 100644
--- a/src/MonoInts/ao_mono_ints.irp.f
+++ b/src/MonoInts/ao_mono_ints.irp.f
@@ -119,3 +119,17 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
   !$OMP END PARALLEL DO
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)]
+  implicit none
+  integer :: i,j,n,l
+  BEGIN_DOC
+ ! array of the mono electronic hamiltonian on the AOs basis
+ ! : sum of the kinetic and nuclear electronic potential 
+  END_DOC
+  do j = 1, ao_num
+   do i = 1, ao_num
+    ao_mono_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + ao_kinetic_integral(i,j)
+   enddo
+  enddo
+END_PROVIDER
+
diff --git a/src/MonoInts/mo_mono_ints.irp.f b/src/MonoInts/mo_mono_ints.irp.f
index 956c6454..e466c9cc 100644
--- a/src/MonoInts/mo_mono_ints.irp.f
+++ b/src/MonoInts/mo_mono_ints.irp.f
@@ -39,8 +39,8 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to
  ! array of the mono electronic hamiltonian on the MOs basis
  ! : sum of the kinetic and nuclear electronic potential 
   END_DOC
-  do i = 1, mo_tot_num
-   do j = 1, mo_tot_num
+  do j = 1, mo_tot_num
+   do i = 1, mo_tot_num
     mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j)
    enddo
   enddo