diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f
index a4166e10..e7480ca2 100644
--- a/src/Determinants/davidson.irp.f
+++ b/src/Determinants/davidson.irp.f
@@ -492,34 +492,42 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
       ! --------------------------------------------------
       
      !$OMP PARALLEL DEFAULT(NONE)                                     &
-        !$OMP PRIVATE(k,i,l,iter2) SHARED(U,W,R,y,iter,lambda,N_st,sze)
+        !$OMP PRIVATE(k,i,l,iter2) &
+        !$OMP SHARED(U,W,R,y,iter,lambda,N_st,sze,to_print, &
+        !$OMP residual_norm,nuclear_repulsion)
       do k=1,N_st
         !$OMP DO
         do i=1,sze
           U(i,k,iter+1) = 0.d0
           W(i,k,iter+1) = 0.d0
+        enddo
+        !$OMP END DO
+         do iter2=1,iter
           do l=1,N_st
-            do iter2=1,iter
+            !$OMP DO
+            do i=1,sze
               U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1)
               W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1)
             enddo
+            !$OMP END DO NOWAIT
           enddo
         enddo
-        !$OMP END DO
-      enddo
-      !$OMP END PARALLEL
       
-      ! Compute residual vector
-      ! -----------------------
+        ! Compute residual vector
+        ! -----------------------
       
-      do k=1,N_st
+        !$OMP DO
         do i=1,sze
           R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1)
         enddo
+        !$OMP END DO
+        !$OMP SINGLE
         residual_norm(k) = u_dot_u(R(1,k),sze)
         to_print(1,k) = lambda(k) + nuclear_repulsion
         to_print(2,k) = residual_norm(k)
+        !$OMP END SINGLE
       enddo
+      !$OMP END PARALLEL
       
       write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))')  iter, to_print(:,1:N_st)
       call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)