diff --git a/plugins/Hartree_Fock/DIIS.irp.f b/plugins/Hartree_Fock/DIIS.irp.f
new file mode 100644
index 00000000..6c5df0f5
--- /dev/null
+++ b/plugins/Hartree_Fock/DIIS.irp.f
@@ -0,0 +1,166 @@
+begin_provider [double precision, FPS_SPF_Matrix_AO, (AO_num_align, AO_num)]
+  implicit none
+  begin_doc 
+!   Commutator FPS - SPF
+  end_doc
+  double precision, allocatable:: scratch(:,:)
+  allocate(                       &
+    scratch(AO_num_align, AO_num) &
+  )
+
+! Compute FP
+
+  call dgemm('N','N',AO_num,AO_num,AO_num,                &
+       1.d0,                                              &
+       Fock_Matrix_AO,Size(Fock_Matrix_AO,1),             &
+       HF_Density_Matrix_AO,Size(HF_Density_Matrix_AO,1), &
+       0.d0,                                              & 
+       scratch,Size(scratch,1))
+
+! Compute FPS
+
+  call dgemm('N','N',AO_num,AO_num,AO_num, &
+       1.d0,                               &
+       scratch,Size(scratch,1),            & 
+       AO_Overlap,Size(AO_Overlap,1),      &
+       0.d0,                               &
+       FPS_SPF_Matrix_AO,Size(FPS_SPF_Matrix_AO,1))
+
+! Compute SP
+
+  call dgemm('N','N',AO_num,AO_num,AO_num,                &
+       1.d0,                                              &
+       AO_Overlap,Size(AO_Overlap,1),                     & 
+       HF_Density_Matrix_AO,Size(HF_Density_Matrix_AO,1), &
+       0.d0,                                              & 
+       scratch,Size(scratch,1))
+
+! Compute FPS - SPF
+
+  call dgemm('N','N',AO_num,AO_num,AO_num,     &
+       -1.d0,                                  & 
+       scratch,Size(scratch,1),                &
+       Fock_Matrix_AO,Size(Fock_Matrix_AO,1),  &
+       1.d0,                                   & 
+       FPS_SPF_Matrix_AO,Size(FPS_SPF_Matrix_AO,1))
+
+end_provider
+
+
+ BEGIN_PROVIDER [ double precision, eigenvalues_Fock_matrix_AO, (AO_num) ]
+&BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_AO, (AO_num_align,AO_num) ]
+
+   BEGIN_DOC
+   ! Eigenvalues and eigenvectors of the Fock matrix over the AO basis
+   END_DOC
+
+   implicit none
+   
+   double precision, allocatable  :: scratch(:,:),work(:),Xt(:,:)
+   integer                        :: lwork,info
+   integer                        :: i,j
+   
+   lwork = 3*AO_num - 1
+   allocate(                       &
+     scratch(AO_num_align,AO_num), &
+     work(lwork),                  &
+     Xt(AO_num,AO_num)             &
+   )
+ 
+! Calculate Xt
+
+  do i=1,AO_num
+    do j=1,AO_num
+      Xt(i,j) = X_Matrix_AO(j,i)
+    enddo
+  enddo
+
+! Calculate Fock matrix in orthogonal basis: F' = Xt.F.X
+
+  call dgemm('N','N',AO_num,AO_num,AO_num,     &
+       1.d0,                                   &
+       Fock_matrix_AO,size(Fock_matrix_AO,1),  &
+       X_Matrix_AO,size(X_Matrix_AO,1),        &
+       0.d0,                                   &
+       eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))       
+
+  call dgemm('N','N',AO_num,AO_num,AO_num,                              &
+       1.d0,                                                            &
+       Xt,size(Xt,1),                                                   &
+       eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1), &
+       0.d0,                                                            &
+       scratch,size(scratch,1))
+     
+! Diagonalize F' to obtain eigenvectors in orthogonal basis C' and eigenvalues
+  
+   call dsyev('V','U',AO_num,       &
+        scratch,size(scratch,1),    &
+        eigenvalues_Fock_matrix_AO, &
+        work,lwork,info)
+ 
+   if(info /= 0) then
+     print *,  irp_here//' failed : ', info
+     stop 1
+   endif
+
+! Back-transform eigenvectors: C =X.C'
+
+  call dgemm('N','N',AO_num,AO_num,AO_num,     &
+       1.d0,                                   &
+       X_matrix_AO,size(X_matrix_AO,1),        &
+       scratch,size(scratch,1),                &
+       0.d0,                                   &
+       eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))       
+   
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num_align,AO_num) ]
+
+  BEGIN_DOC
+!   Matrix X = S^{-1/2} obtained by SVD
+  END_DOC
+
+  implicit none
+  
+  integer                         :: LDA, LDC
+  double precision, allocatable   :: U(:,:),Vt(:,:), D(:)
+  integer                         :: info, i, j, k
+ 
+  LDA = size(AO_overlap,1)
+  LDC = size(AO_overlap,1)
+ 
+  allocate(         &
+    U(LDC,AO_num),  &
+    Vt(LDA,AO_num), &
+    D(AO_num))
+
+  call svd(            &
+       AO_overlap,LDA, &
+       U,LDC,          &
+       D,              &
+       Vt,LDA,         &
+       AO_num,AO_num)
+
+  do i=1,AO_num
+    if(abs(D(i)) < threshold_overlap_AO_eigenvalues) then
+      D(i) = 0.d0
+    else
+      D(i) = 1.d0/sqrt(D(i))
+    endif
+    do j=1,AO_num
+      X_matrix_AO(j,i) = 0.d0
+    enddo
+  enddo
+
+  do k=1,AO_num
+    if(D(k) /= 0.d0) then
+      do j=1,AO_num
+        do i=1,MO_tot_num
+          X_matrix_AO(i,j) = X_matrix_AO(i,j) + U(i,k)*D(k)*Vt(k,j)
+        enddo
+      enddo
+    endif
+  enddo
+  
+
+END_PROVIDER
diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg
index 2fa29cf0..f8923e57 100644
--- a/plugins/Hartree_Fock/EZFIO.cfg
+++ b/plugins/Hartree_Fock/EZFIO.cfg
@@ -1,3 +1,21 @@
+[threshold_overlap_ao_eigenvalues]
+type: Threshold
+doc: Threshold on the magnitude of the smallest eigenvalues of the overlap matrix in the AO basis
+interface: ezfio,provider,ocaml
+default: 1.e-6
+
+[max_dim_diis]
+type: integer
+doc: Maximum size of the DIIS extrapolation procedure
+interface: ezfio,provider,ocaml
+default: 15
+
+[threshold_diis]
+type: Threshold
+doc: Threshold on the convergence of the DIIS error vector during a Hartree-Fock calculation
+interface: ezfio,provider,ocaml
+default: 1.e-5
+
 [thresh_scf]
 type: Threshold
 doc: Threshold on the convergence of the Hartree Fock energy
@@ -8,7 +26,7 @@ default: 1.e-10
 type: Strictly_positive_int
 doc: Maximum number of SCF iterations
 interface: ezfio,provider,ocaml
-default: 200
+default: 128
 
 [level_shift]
 type: Positive_float
@@ -16,6 +34,12 @@ doc: Energy shift on the virtual MOs to improve SCF convergence
 interface: ezfio,provider,ocaml
 default: 0.5
 
+[scf_algorithm]
+type: character*(32)
+doc: Type of SCF algorithm used. Possible choices are [ damp | DIIS]
+interface: ezfio,provider,ocaml
+default: damp
+
 [mo_guess_type]
 type: MO_guess
 doc: Initial MO guess. Can be [ Huckel | HCore ]
diff --git a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
new file mode 100644
index 00000000..989ae387
--- /dev/null
+++ b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
@@ -0,0 +1,218 @@
+subroutine Roothaan_Hall_SCF
+
+BEGIN_DOC
+! Roothaan-Hall algorithm for SCF Hartree-Fock calculation
+END_DOC
+
+  implicit none
+
+  double precision               :: energy_SCF,energy_SCF_previous,Delta_energy_SCF,max_error_DIIS
+  double precision, allocatable  :: Fock_matrix_DIIS(:,:,:),error_matrix_DIIS(:,:,:)
+
+  integer                        :: iteration_SCF,dim_DIIS,index_dim_DIIS
+ 
+  integer                        :: i,j
+ 
+  allocate(                                       &
+    Fock_matrix_DIIS(AO_num,AO_num,max_dim_DIIS), &
+    error_matrix_DIIS(AO_num,AO_num,max_dim_DIIS) &
+  )
+
+  call write_time(output_hartree_fock)
+
+  write(output_hartree_fock,'(A4, 1X, A16, 1X, A16, 1X, A16)')  &
+    '====','================','================','================'
+  write(output_hartree_fock,'(A4, 1X, A16, 1X, A16, 1X, A16)')  &
+    '  N ', 'Energy  ', 'Energy diff  ', 'DIIS error  '
+  write(output_hartree_fock,'(A4, 1X, A16, 1X, A16, 1X, A16)')  &
+    '====','================','================','================'
+
+! Initialize energies and density matrices
+
+  energy_SCF_previous = HF_energy
+  Delta_energy_SCF = 0.d0
+  iteration_SCF = 0
+  dim_DIIS = 0
+  max_error_DIIS = 1.d0
+
+!
+! Start of main SCF loop
+!
+  do while((max_error_DIIS > threshold_DIIS) .and. (iteration_SCF < n_it_SCF_max))
+
+! Increment cycle number
+
+    iteration_SCF += 1
+
+! Current size of the DIIS space
+
+    dim_DIIS = min(dim_DIIS+1,max_dim_DIIS)
+ 
+! Store Fock and error matrices at each iteration
+    
+    do j=1,AO_num
+      do i=1,AO_num
+        index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
+        Fock_matrix_DIIS(i,j,index_dim_DIIS) = Fock_matrix_AO(i,j)
+        error_matrix_DIIS(i,j,index_dim_DIIS) = FPS_SPF_Matrix_AO(i,j)
+      enddo
+    enddo
+    
+! Compute the extrapolated Fock matrix
+
+  call extrapolate_Fock_matrix(         &
+    error_matrix_DIIS,Fock_matrix_DIIS, &
+    iteration_SCF,dim_DIIS              &
+  ) 
+
+    touch Fock_matrix_AO
+
+    MO_coef = eigenvectors_Fock_matrix_AO
+
+! This algorithm still have an issue with linear dependencies
+!   do i=1,AO_num
+!     write(*,*) i,eigenvalues_Fock_matrix_AO(i)
+!   enddo
+
+    touch MO_coef
+
+! Calculate error vectors
+
+    max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_AO))
+
+!   SCF energy
+
+    energy_SCF = HF_energy
+    Delta_Energy_SCF = energy_SCF - energy_SCF_previous
+    energy_SCF_previous = energy_SCF
+
+! Print results at the end of each iteration
+
+    write(output_hartree_fock,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10)')  &
+      iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS
+
+  enddo
+
+!
+! End of Main SCF loop
+!
+
+  write(output_hartree_fock,'(A4, 1X, A16, 1X, A16, 1X, A16)') &
+    '====','================','================','================'
+  write(output_hartree_fock,*)
+  
+  if(.not.no_oa_or_av_opt)then
+   call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
+  endif
+
+  call write_double(output_hartree_fock, Energy_SCF, 'Hartree-Fock energy')
+  call ezfio_set_hartree_fock_energy(Energy_SCF)
+
+  call write_time(output_hartree_fock)
+
+end
+
+subroutine extrapolate_Fock_matrix(   &
+  error_matrix_DIIS,Fock_matrix_DIIS, &
+  iteration_SCF,dim_DIIS              &
+) 
+
+BEGIN_DOC
+! Compute the extrapolated Fock matrix using the DIIS procedure
+END_DOC
+
+  implicit none
+
+  double precision,intent(in)   :: Fock_matrix_DIIS(AO_num,AO_num,*),error_matrix_DIIS(AO_num,AO_num,*)
+  integer,intent(in)            :: iteration_SCF
+  integer,intent(inout)         :: dim_DIIS
+
+  double precision,allocatable  :: B_matrix_DIIS(:,:),X_vector_DIIS(:)
+
+  double precision,allocatable  :: scratch(:,:)
+  integer                       :: i,j,k,i_DIIS,j_DIIS
+
+  allocate(                               &
+    B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), &
+    X_vector_DIIS(dim_DIIS+1),            &
+    scratch(AO_num,AO_num)                &
+  )
+
+! Compute the matrices B and X
+  do j=1,dim_DIIS
+    do i=1,dim_DIIS
+
+      j_DIIS = mod(iteration_SCF-j,max_dim_DIIS)+1
+      i_DIIS = mod(iteration_SCF-i,max_dim_DIIS)+1
+
+! Compute product of two errors vectors
+
+      call dgemm('N','N',AO_num,AO_num,AO_num,                      &
+           1.d0,                                                    &
+           error_matrix_DIIS(1,1,i_DIIS),size(error_matrix_DIIS,1), &
+           error_matrix_DIIS(1,1,j_DIIS),size(error_matrix_DIIS,1), &
+           0.d0,                                                    &
+           scratch,size(scratch,1))                                                  
+
+! Compute Trace
+
+      B_matrix_DIIS(i,j) = 0.d0
+      do k=1,AO_num
+        B_matrix_DIIS(i,j) += scratch(k,k)
+      enddo
+    enddo
+  enddo
+
+! Pad B matrix and build the X matrix
+
+  do i=1,dim_DIIS
+    B_matrix_DIIS(i,dim_DIIS+1) = -1.d0
+    B_matrix_DIIS(dim_DIIS+1,i) = -1.d0
+    X_vector_DIIS(i) = 0.d0
+  enddo
+  B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1) = 0.d0
+  X_vector_DIIS(dim_DIIS+1) = -1.d0
+
+! Solve the linear system C = B.X
+
+  integer              :: info
+  integer,allocatable  :: ipiv(:)
+
+  allocate(          &
+    ipiv(dim_DIIS+1) &
+  )
+
+  call dsysv('U',dim_DIIS+1,1,           &
+    B_matrix_DIIS,size(B_matrix_DIIS,1), &
+    ipiv,                                &
+    X_vector_DIIS,size(X_vector_DIIS,1), &
+    scratch,size(scratch),               &
+    info                                 &
+  )
+
+ if(info == 0) then
+   
+! Compute extrapolated Fock matrix
+
+    Fock_matrix_AO(:,:) = 0.d0
+
+    do k=1,dim_DIIS
+      do j=1,AO_num
+        do i=1,AO_num
+          Fock_matrix_AO(i,j) += X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1)
+        enddo
+      enddo
+    enddo
+
+  else
+    write(*,*) 'Re-initialize DIIS!!'
+    dim_DIIS = 0
+  endif
+
+! do i=1,AO_num
+!   do j=1,AO_num
+!     write(*,*) Fock_matrix_AO(i,j)
+!   enddo
+! enddo
+
+end
diff --git a/plugins/Hartree_Fock/SCF.irp.f b/plugins/Hartree_Fock/SCF.irp.f
index dead61ee..715bea86 100644
--- a/plugins/Hartree_Fock/SCF.irp.f
+++ b/plugins/Hartree_Fock/SCF.irp.f
@@ -13,7 +13,7 @@ end
 subroutine create_guess
   implicit none
   BEGIN_DOC
-! Create an MO guess if no MOs are present in the EZFIO directory
+!   Create a MO guess if no MOs are present in the EZFIO directory
   END_DOC
   logical                        :: exists
   PROVIDE ezfio_filename
@@ -34,21 +34,34 @@ subroutine create_guess
   endif
 end
 
+ao_to_mo
 
 subroutine run
 
+  BEGIN_DOC
+!   Run SCF calculation
+  END_DOC
+
   use bitmasks
   implicit none
-  BEGIN_DOC
-! Run SCF calculation
-  END_DOC
+
   double precision               :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem
-  double precision               :: E0
+  double precision               :: EHF
   integer                        :: i_it, i, j, k
    
-  E0 = HF_energy 
+  EHF = HF_energy 
 
   mo_label = "Canonical"
-  call damping_SCF
+
+! Choose SCF algorithm
+
+  if(scf_algorithm == 'damp') then
+    call damping_SCF
+  else if(scf_algorithm == 'DIIS') then
+    call Roothaan_Hall_SCF
+  else
+    write(*,*) 'Unrecognized SCF algorithm: '//scf_algorithm
+    stop 1
+  endif
   
 end