diff --git a/src/nloptics/el_density_vbands.F b/src/nloptics/el_density_vbands.F
index 896f4b373..b786db98b 100644
--- a/src/nloptics/el_density_vbands.F
+++ b/src/nloptics/el_density_vbands.F
@@ -47,23 +47,22 @@ subroutine el_density_vbands(en,Xk,rho,VAL_BANDS)
  integer     :: i1,ik,is,rho_syms,ir,ik_mem
  real(SP)    :: rho_no_sym(fft_size)
  !
-!$OMP WORKSHARE  
- rho=rZERO
  rho_no_sym=rZERO
-!$OMP END WORKSHARE 
  !
  do ik=1,Xk%nibz
    !  
    if (.not.PAR_IND_Xk_ibz%element_1D(ik)) cycle
    ik_mem=PAR_Xk_ibz_index(ik)
    !
+   !$omp parallel do, default(shared), private(is,ir), &
+   !$omp              reduction(+:rho_no_sym), collapse(2)
    do is=1,n_sp_pol
-!$OMP WORKSHARE  
-     forall(ir=1:fft_size)
-       rho_no_sym(ir)=rho_no_sym(ir)+real(spin_occ,SP)*Xk%weights(ik)*sum(abs(VAL_BANDS(ir,:en%nbf(is),ik_mem,is))**2._SP)
-     end forall
-!$OMP END WORKSHARE  
+     do ir=1,fft_size
+       rho_no_sym(ir)=rho_no_sym(ir)+real(spin_occ,SP)*Xk%weights(ik)*&
+&                         sum(abs(VAL_BANDS(ir,:en%nbf(is),ik_mem,is))**2._SP)
+     enddo
    enddo
+   !$omp end parallel do
    !
  enddo
  !
@@ -71,18 +70,22 @@ subroutine el_density_vbands(en,Xk,rho,VAL_BANDS)
  !
  ! Simmetrization
  !
-!$OMP WORKSHARE  
  rho_syms=nsym/(i_time_rev+1)
-!$OMP END WORKSHARE  
- !
-!$OMP WORKSHARE  
- forall(ir=1:fft_size)
+ rho(:)=rZERO
+ ! 
+!$omp parallel default(shared), private(ir)
+!$omp do, reduction(+:rho)
+ do ir=1,fft_size
    rho(ir)=rho(ir)+sum(rho_no_sym(fft_rot_r(ir,1:rho_syms)))/real(nsym,SP)
- end forall
-!$OMP END WORKSHARE  
+ enddo
+!$omp end do
  !
-!$OMP WORKSHARE  
- rho=real(1._SP+i_time_rev,SP)*rho
-!$OMP END WORKSHARE  
+!$omp do
+ do ir=1,fft_size
+   rho(ir)=real(1._SP+i_time_rev,SP)*rho(ir)
+ enddo
+!$omp end do
+!$omp end parallel
  !
 end subroutine
+