diff --git a/XClib/qe_drivers_mgga.f90 b/XClib/qe_drivers_mgga.f90
index ae7022af1..cb63e559c 100644
--- a/XClib/qe_drivers_mgga.f90
+++ b/XClib/qe_drivers_mgga.f90
@@ -79,6 +79,9 @@ SUBROUTINE tau_xc( length, rho, grho2, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c
     IF ( (arho<=rho_threshold_mgga).OR.(grho2(k)<=grho2_threshold_mgga).OR. &
          (ABS(tau(k))<=rho_threshold_mgga) ) CYCLE
     !
+    ! ...libxc-like threshold management
+    !grho2(k) = MIN( grho2(k), (8.d0*rho(k)*tau(k))**2 )
+    !
     SELECT CASE( imeta )
     CASE( 1 )
        !
@@ -147,14 +150,15 @@ SUBROUTINE tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c,
   ex=0.0_DP ; v1x=0.0_DP ; v2x=0.0_DP ; v3x=0.0_DP
   ec=0.0_DP ; v1c=0.0_DP ; v2c=0.0_DP ; v3c=0.0_DP
   !
-  ! FIXME: for SCAN, this will be calculated later
-  !
   DO k = 1, length
      !
      rh   = rho(k,1) + rho(k,2)
      atau = tau(k,1) + tau(k,2)             ! KE-density in Hartree
-     grho2(1) = SUM( grho(:,k,1)**2 )
-     grho2(2) = SUM( grho(:,k,2)**2 )
+     ! ...libxc-like threshold management
+     !grho2(1) = MIN( SUM(grho(:,k,1)**2), (8.d0*rho(k,1)*tau(k,1))**2 )
+     !grho2(2) = MIN( SUM(grho(:,k,2)**2), (8.d0*rho(k,2)*tau(k,2))**2 )
+     grho2(1) = SUM(grho(:,k,1)**2) 
+     grho2(2) = SUM(grho(:,k,2)**2)
      ggrho2 = ( grho2(1) + grho2(2) ) * 4.0_DP
      !
      IF ( (rh <= rho_threshold_mgga).OR.(ggrho2 <= grho2_threshold_mgga).OR.&
diff --git a/XClib/qe_funct_mgga.f90 b/XClib/qe_funct_mgga.f90
index 523bed65f..004fccc80 100644
--- a/XClib/qe_funct_mgga.f90
+++ b/XClib/qe_funct_mgga.f90
@@ -267,9 +267,7 @@ SUBROUTINE metac( rho, grho2, tau, ec, v1c, v2c, v3c )                    !<GPU:
   cf3 = cf3*rho
   v2c = cf2*d2rev + cf3*2.0_DP/grho2
   v3c = cf2*d3rev - cf3/tau
-  !
   ec = rho*ec_rev*(1.0_DP+dd*ec_rev*z2*z)  !-rho*ec_unif(1)
-  !v1c = v1c - vc_unif(1)
   !
   RETURN
   !
@@ -865,8 +863,8 @@ END SUBROUTINE metac_spin
 !
 !           input:  - rho
 !                   - grho2=|\nabla rho|^2
-!                   - tau = the kinetic energy density
-!                           It is defined as summ_i( |nabla phi_i|**2 )
+!                   - tau = the input kinetic energy density
+!                           It is defined as 0.5summ_i( |nabla phi_i|**2 )
 !
 !           definition:  E_x = \int ex dr
 !
@@ -923,18 +921,20 @@ SUBROUTINE m06lxc( rho, grho2, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c )
   !
   taua = tau * two * 0.5_dp ! Taua, which is Tau_sigma is half Tau
   taub = taua               ! Tau is defined as summ_i( |nabla phi_i|**2 )
-                              ! in the M06L routine
+                              ! in the following M06L routines
   !
   CALL m06lx( rhoa, grho2a, taua, ex, v1x, v2x, v3x )   !<GPU:m06lx=>m06lx_d>
   !
   ex = two * ex  ! Add the two components up + dw
   !
   v2x = 0.5_dp * v2x
+  v3x = 2.0_dp * v3x
   !
   CALL m06lc( rhoa, rhob, grho2a, grho2b, taua, taub, ec, v1c, v2c, v3c, &   !<GPU:m06lc=>m06lc_d>
               v1cb, v2cb, v3cb )
   !
   v2c = 0.5_dp * v2c
+  v3c = 2.0_dp * v3c
   !
 END SUBROUTINE m06lxc
 !
@@ -963,9 +963,13 @@ SUBROUTINE m06lxc_spin( rhoup, rhodw, grhoup2, grhodw2, tauup, taudw,      &
   CALL m06lx( rhodw, grhodw2, taub, exdw, v1xdw, v2xdw, v3xdw )   !<GPU:m06lx=>m06lx_d>
   !
   ex = exup + exdw
+  v3xup = 2.0_dp * v3xup
+  v3xdw = 2.0_dp * v3xdw
   !
   CALL m06lc( rhoup, rhodw, grhoup2, grhodw2, taua, taub, &       !<GPU:m06lc=>m06lc_d>
               ec, v1cup, v2cup, v3cup, v1cdw, v2cdw, v3cdw )
+  v3cup = 2.0_dp * v3cup
+  v3cdw = 2.0_dp * v3cdw
   !
 END SUBROUTINE m06lxc_spin
 ! !
diff --git a/XClib/xc_wrapper_gga.f90 b/XClib/xc_wrapper_gga.f90
index 1f518fa7d..4d45187c1 100644
--- a/XClib/xc_wrapper_gga.f90
+++ b/XClib/xc_wrapper_gga.f90
@@ -174,19 +174,23 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
       ENDDO
     ELSE
       DO k = 1, length
-        sgn(:) = 1.d0
-        !IF (rho_lxc(2*k-1)<rho_threshold_gga .OR. SQRT(ABS(sigma(3*k-2)))<grho_threshold_gga) sgn(1)=0.d0
-        !IF (rho_lxc(2*k)  <rho_threshold_gga .OR. SQRT(ABS(sigma(3*k)))  <grho_threshold_gga) sgn(2)=0.d0
-        IF (rho_lxc(2*k-1)<rho_threshold_gga) sgn(1)=0.d0
-        !.OR. SQRT(ABS(sigma(3*k-2)))<grho_threshold_gga) sgn(1)=0.d0
-        IF (rho_lxc(2*k)  <rho_threshold_gga) sgn(2)=0.d0
-        !.OR. SQRT(ABS(sigma(3*k)))  <grho_threshold_gga) sgn(2)=0.d0
-        ec(k) = ec_lxc(k) * (rho_lxc(2*k-1)*sgn(1)+rho_lxc(2*k)*sgn(2))
-        v1c(k,1) = vc_rho(2*k-1) * sgn(1)
-        v1c(k,2) = vc_rho(2*k) * sgn(2)
-        v2c(k,1) = vc_sigma(3*k-2)*2.d0 * sgn(1)
-        v2c_ud(k)= vc_sigma(3*k-1) * sgn(1)*sgn(2)
-        v2c(k,2) = vc_sigma(3*k)*2.d0 * sgn(2)
+        rho_up = rho_lxc(2*k-1)
+        rho_dw = rho_lxc(2*k)
+        grho2_up = sigma(3*k-2)
+        grho2_dw = sigma(3*k)
+        xlda_up = 1.0_DP ; xlda_dw = 1.0_DP
+        xgga_up = 1.0_DP ; xgga_dw = 1.0_DP
+        IF ( rho_up+rho_dw <= rho_threshold_gga ) CYCLE
+        IF ( rho_up <= rho_threshold_lda ) xlda_up = 0.0_DP
+        IF ( rho_dw <= rho_threshold_lda ) xlda_dw = 0.0_DP
+        IF ( rho_up<=small .OR. SQRT(ABS(grho2_up))<=small ) xgga_up = 0.0_DP
+        IF ( rho_dw<=small .OR. SQRT(ABS(grho2_dw))<=small ) xgga_dw = 0.0_DP
+        ec(k) = ec_lxc(k) * (rho_up*xlda_up+rho_dw*xlda_dw)
+        v1c(k,1) = vc_rho(2*k-1)*xlda_up
+        v1c(k,2) = vc_rho(2*k)*xlda_dw
+        v2c(k,1) = vc_sigma(3*k-2)*2.d0*xgga_up
+        v2c_ud(k)= vc_sigma(3*k-1)*xgga_up*xgga_dw
+        v2c(k,2) = vc_sigma(3*k)*2.d0*xgga_dw
       ENDDO
     ENDIF
     !  
diff --git a/XClib/xc_wrapper_mgga.f90 b/XClib/xc_wrapper_mgga.f90
index c2839acab..7181db17c 100644
--- a/XClib/xc_wrapper_mgga.f90
+++ b/XClib/xc_wrapper_mgga.f90
@@ -19,9 +19,9 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
 #endif 
   !
   USE kind_l,               ONLY: DP
-  USE dft_setting_params,   ONLY: imeta, imetac, is_libxc, rho_threshold_mgga,  &
-                            grho2_threshold_mgga, tau_threshold_mgga, scan_exx, &
-                            exx_started, exx_fraction
+  USE dft_setting_params,   ONLY: imetac, is_libxc, rho_threshold_mgga,     &
+                                  grho2_threshold_mgga, tau_threshold_mgga, &
+                                  scan_exx, exx_started, exx_fraction
   USE qe_drivers_mgga
   !
   IMPLICIT NONE
@@ -59,16 +59,18 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
   !
   INTEGER :: is
   REAL(DP), ALLOCATABLE :: grho2(:,:)
+  REAL(DP), PARAMETER :: small = 1.E-10_DP
   !
 #if defined(__LIBXC)
   REAL(DP), ALLOCATABLE :: rho_lxc(:), sigma(:), tau_lxc(:)
   REAL(DP), ALLOCATABLE :: ex_lxc(:), ec_lxc(:)
   REAL(DP), ALLOCATABLE :: vx_rho(:), vx_sigma(:), vx_tau(:)
   REAL(DP), ALLOCATABLE :: vc_rho(:), vc_sigma(:), vc_tau(:)
-  REAL(DP), ALLOCATABLE :: lapl_rho(:), vlapl_rho(:) ! not used in TPSS
+  REAL(DP), ALLOCATABLE :: lapl_rho(:), vlapl_rho(:) ! not used in QE
   !
-  INTEGER :: k, ipol, pol_unpol, eflag
-  LOGICAL :: POLARIZED
+  INTEGER  :: k, ipol, pol_unpol, eflag
+  LOGICAL  :: POLARIZED
+  REAL(DP) :: rh, ggrho2, atau
 #if (XC_MAJOR_VERSION > 4)
   INTEGER(8) :: lengthxc
 #else
@@ -103,6 +105,7 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
       sigma(k) = MAX( grho(1,k,1)**2 + grho(2,k,1)**2 + grho(3,k,1)**2, &
                       grho2_threshold_mgga )
       tau_lxc(k) = MAX( tau(k,1), tau_threshold_mgga )
+      vlapl_rho(k) = 0._DP
     ENDDO
     !
   ELSE
@@ -110,16 +113,16 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
     DO k = 1, length
        rho_lxc(2*k-1) = ABS( rho(k,1) )
        rho_lxc(2*k)   = ABS( rho(k,2) )
-       !
        sigma(3*k-2) = MAX( grho(1,k,1)**2 + grho(2,k,1)**2 + grho(3,k,1)**2, &
                            grho2_threshold_mgga )
-       sigma(3*k-1) = grho(1,k,1) * grho(1,k,2) + grho(2,k,1) * grho(2,k,2) + &
+       sigma(3*k-1) = grho(1,k,1) * grho(1,k,2) + grho(2,k,1) * grho(2,k,2) +&
                       grho(3,k,1) * grho(3,k,2)
        sigma(3*k)   = MAX( grho(1,k,2)**2 + grho(2,k,2)**2 + grho(3,k,2)**2, &
                            grho2_threshold_mgga )
-       !
-       tau_lxc(2*k-1) = MAX( tau(k,1), tau_threshold_mgga )
-       tau_lxc(2*k)   = MAX( tau(k,2), tau_threshold_mgga )
+       tau_lxc(2*k-1) = MAX( tau(k,1), small )
+       tau_lxc(2*k)   = MAX( tau(k,2), small )
+       vlapl_rho(2*k-1) = 0._DP
+       vlapl_rho(2*k)   = 0._DP
     ENDDO
     !
   ENDIF
@@ -132,19 +135,17 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
        grho2(:,is) = grho(1,:,is)**2 + grho(2,:,is)**2 + grho(3,:,is)**2
     ENDDO
     !
-    IF (ns == 1) THEN
-       CALL tau_xc( length, rho(:,1), grho2(:,1), tau(:,1), ex, ec, v1x(:,1), &
-                    v2x(:,1), v3x(:,1), v1c(:,1), v2c(1,:,1), v3c(:,1) )
-    ELSEIF (ns == 2) THEN
-       CALL tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, &
-                         v2c, v3c )
-    ENDIF
+    IF (ns == 1) CALL tau_xc( length, rho(:,1), grho2(:,1), tau(:,1), ex, ec,    &
+                              v1x(:,1), v2x(:,1), v3x(:,1), v1c(:,1), v2c(1,:,1),&
+                              v3c(:,1) )
+    IF (ns == 2) CALL tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x,&
+                                   v1c, v2c, v3c )
     !
     DEALLOCATE( grho2 )
     !
   ENDIF
   !
-  ! META EXCHANGE
+  ! ... META EXCHANGE
   !
   IF ( is_libxc(5) ) THEN
     CALL xc_f03_func_set_dens_threshold( xc_func(5), rho_threshold_mgga )
@@ -160,6 +161,9 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
     !
     IF (.NOT. POLARIZED) THEN
       DO k = 1, length
+        IF ( ABS(rho_lxc(k))<=rho_threshold_mgga .OR. &
+             sigma(k)<=grho2_threshold_mgga      .OR. &
+             ABS(tau_lxc(k))<=rho_threshold_mgga ) CYCLE
         ex(k) = ex_lxc(k) * rho_lxc(k)
         v1x(k,1) = vx_rho(k)
         v2x(k,1) = vx_sigma(k) * 2.0_DP
@@ -167,13 +171,22 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
       ENDDO
     ELSE
       DO k = 1, length
+        IF (rho_lxc(2*k-1)+rho_lxc(2*k) <= rho_threshold_mgga) CYCLE
         ex(k) = ex_lxc(k) * (rho_lxc(2*k-1)+rho_lxc(2*k))
-        v1x(k,1) = vx_rho(2*k-1)
-        v1x(k,2) = vx_rho(2*k)
-        v2x(k,1) = vx_sigma(3*k-2)*2.d0
-        v2x(k,2) = vx_sigma(3*k)*2.d0
-        v3x(k,1) = vx_tau(2*k-1)
-        v3x(k,2) = vx_tau(2*k)
+        IF ( ABS(rho_lxc(2*k-1))>rho_threshold_mgga .AND. &
+             sigma(3*k-2)>grho2_threshold_mgga      .AND. &
+             ABS(tau_lxc(2*k-1))>tau_threshold_mgga ) THEN
+          v1x(k,1) = vx_rho(2*k-1)
+          v2x(k,1) = vx_sigma(3*k-2)*2.d0
+          v3x(k,1) = vx_tau(2*k-1)
+        ENDIF
+        IF ( ABS(rho_lxc(2*k))>rho_threshold_mgga .AND. &
+             sigma(3*k)>grho2_threshold_mgga      .AND. &
+             ABS(tau_lxc(2*k))>tau_threshold_mgga ) THEN
+          v1x(k,2) = vx_rho(2*k)
+          v2x(k,2) = vx_sigma(3*k)*2.d0
+          v3x(k,2) = vx_tau(2*k)                                   !m06l: *2 libxc o /2 qe????
+        ENDIF                                                      !-chiarisci il problema del fattore 2
       ENDDO
     ENDIF
     !
@@ -189,7 +202,7 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
     !
   ENDIF
   !
-  ! META CORRELATION
+  ! ... META CORRELATION
   !
   IF ( is_libxc(6) ) THEN
     !
@@ -206,13 +219,22 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
     !
     IF (.NOT. POLARIZED) THEN
        DO k = 1, length
-          ec(k) = ec_lxc(k) * rho_lxc(k) 
-          v1c(k,1) = vc_rho(k)
-          v2c(1,k,1) = vc_sigma(k) * 2.0_DP
-          v3c(k,1) = vc_tau(k)
+         IF ( ABS(rho_lxc(k))<=rho_threshold_mgga    .OR. &
+                     sigma(k)<=grho2_threshold_mgga  .OR. &
+              ABS(tau_lxc(k))<=rho_threshold_mgga  ) CYCLE
+         ec(k) = ec_lxc(k) * rho_lxc(k) 
+         v1c(k,1) = vc_rho(k)
+         v2c(1,k,1) = vc_sigma(k) * 2.0_DP
+         v3c(k,1) = vc_tau(k)
        ENDDO
     ELSE
        DO k = 1, length
+          rh   = rho_lxc(2*k-1) + rho_lxc(2*k)
+          atau = ABS(tau_lxc(2*k-1) + tau_lxc(2*k))
+          ggrho2 = (sigma(3*k-2) + sigma(3*k))*4.0_DP
+          IF ( rh <= rho_threshold_mgga   .OR. &
+           ggrho2 <= grho2_threshold_mgga .OR. &
+             atau <= tau_threshold_mgga  ) CYCLE
           ec(k) = ec_lxc(k) * (rho_lxc(2*k-1)+rho_lxc(2*k))
           v1c(k,1) = vc_rho(2*k-1)
           v1c(k,2) = vc_rho(2*k)
diff --git a/XClib/xclib_test.f90 b/XClib/xclib_test.f90
index fcb746d3a..7fb11eceb 100644
--- a/XClib/xclib_test.f90
+++ b/XClib/xclib_test.f90
@@ -433,7 +433,7 @@ PROGRAM xclib_test
     !
     IF ( TRIM(dft)=='xxxx' .OR.TRIM(dft)=='NONE'   .OR. &
          TRIM(dft)=='META'   .OR. & !TRIM(dft)=='TB09' .OR. &
-         TRIM(dft)=='SCA0' .OR.TRIM(dft)=='TPSS'   .OR. &
+         TRIM(dft)=='SCA0' .OR. & !TRIM(dft)=='TPSS'   .OR. &
          TRIM(dft)=='SCAN0'.OR.TRIM(dft)=='PZ+META'.OR. &
          TRIM(dft)=='PBE+META' ) THEN
          CALL print_test_status( skipped )
@@ -1469,10 +1469,10 @@ PROGRAM xclib_test
     IF (what(1:1)=='V') thr = diff_thr_vmgga
   ENDIF
   !
-  !IF (mype==root .AND. test/='gen-benchmark') THEN
-  !  WRITE(stdout,*) " "
-  !  IF ( POLARIZED .AND. what(1:1)/='E' ) WRITE(stdout,*) TRIM(what)
-  !ENDIF
+  IF (mype==root .AND. TRIM(test)=='exe-benchmark') THEN
+    WRITE(stdout,*) " "
+    IF ( POLARIZED .AND. what(1:1)/='E' ) WRITE(stdout,*) TRIM(what)
+  ENDIF
   !
   xc_aver=0._DP
   !
@@ -1545,11 +1545,6 @@ PROGRAM xclib_test
   IF (LDA .AND. .NOT.GGA) thr = diff_thr_dmuxc
   IF ( GGA ) thr = diff_thr_dv
   !
-  !IF (test=='exe-benchmark') THEN
-  !  WRITE(stdout,*) " "
-  !  WRITE(stdout,*) what
-  !ENDIF
-  !
   dxc_aver=0._DP
   !
   dxc_aver(1,1) = SUM(dxc_qe(1:nnr,1,1))/DBLE(npoints)