diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f
index a63b46b6..3b9eaeb4 100644
--- a/src/scf_utils/roothaan_hall_scf.irp.f
+++ b/src/scf_utils/roothaan_hall_scf.irp.f
@@ -202,7 +202,7 @@ END_DOC
   double precision,allocatable  :: C_vector_DIIS(:)
 
   double precision,allocatable  :: scratch(:,:)
-  integer                       :: i,j,k,i_DIIS,j_DIIS
+  integer                       :: i,j,k,l,i_DIIS,j_DIIS
   double precision :: rcond, ferr, berr
   integer, allocatable :: iwork(:)
   integer :: lwork
@@ -222,38 +222,21 @@ END_DOC
   B_matrix_DIIS(:,:) = 0.d0
   do j=1,dim_DIIS
     j_DIIS = min(dim_DIIS,mod(iteration_SCF-j,max_dim_DIIS)+1)
-    do i=1,dim_DIIS
 
+    do i=1,dim_DIIS
       i_DIIS = min(dim_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
-
-      do k=1,ao_num
-        B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + scratch(k,k)
+      do l=1,ao_num
+        do k=1,ao_num
+          B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + &
+            error_matrix_DIIS(k,l,i_DIIS) * error_matrix_DIIS(k,l,j_DIIS)
+        enddo
       enddo
+
     enddo
   enddo
 
-! TODO : Could be simplified (to be checked)
-!
-!  call dgemm('T','N', &
-!    min(dim_DIIS,iteration_SCF), min(dim_DIIS,iteration_SCF), &
-!    ao_num*ao_num, &
-!    1.d0,                                                    &
-!    error_matrix_DIIS,size(error_matrix_DIIS,1)*size(error_matrix_DIIS,2), &
-!    error_matrix_DIIS,size(error_matrix_DIIS,1)*size(error_matrix_DIIS,2), &
-!    0.d0,                                                    &
-!    B_matrix_DIIS,size(B_matrix_DIIS,1))
-
 ! Pad B matrix and build the X matrix
 
   C_vector_DIIS(:) = 0.d0