diff --git a/CPV/cg_sub.f90 b/CPV/cg_sub.f90
index c3141f8c2..18325f8eb 100644
--- a/CPV/cg_sub.f90
+++ b/CPV/cg_sub.f90
@@ -35,12 +35,11 @@
       use reciprocal_vectors, only: gstart
       use cvan, only: nvb, ish
       use ions_base, only: na, nat, pmass, nax, nsp, rcmax
-      use grid_dimensions, only: nnr => nrxx, nr1, nr2, nr3
+      use grid_dimensions, only: nrxx, nr1, nr2, nr3
       use cell_base, only: ainv, a1, a2, a3
       use cell_base, only: omega, alat
       use cell_base, only: h, hold, deth, wmass, tpiba2
       use smooth_grid_dimensions, only: nrxxs
-      use smallbox_grid_dimensions, only: nnrb => nnrbx, nr1b, nr2b, nr3b
       use local_pseudo, only: vps, rhops
       use io_global,                ONLY : io_global_start, stdout, ionode, ionode_id
       use mp_global,                ONLY : intra_image_comm, np_ortho, me_ortho, ortho_comm, me_image
@@ -80,11 +79,11 @@
       real(dp) :: becdr(nhsa,nspin*nlax,3)
       integer irb(3,nat)
       complex(dp) :: eigrb(ngb,nat)
-      real(dp) :: rhor(nnr,nspin)
-      real(dp) :: vpot(nnr,nspin)
+      real(dp) :: rhor(nrxx,nspin)
+      real(dp) :: vpot(nrxx,nspin)
       complex(dp) :: rhog(ngm,nspin)
       real(dp) :: rhos(nrxxs,nspin)
-      real(dp) :: rhoc(nnr)
+      real(dp) :: rhoc(nrxx)
       complex(dp) :: ei1(-nr1:nr1,nat)
       complex(dp) :: ei2(-nr2:nr2,nat)
       complex(dp) :: ei3(-nr3:nr3,nat)
diff --git a/CPV/chargedensity.f90 b/CPV/chargedensity.f90
index 30bceeb8a..21a641e69 100644
--- a/CPV/chargedensity.f90
+++ b/CPV/chargedensity.f90
@@ -768,11 +768,11 @@ SUBROUTINE drhov(irb,eigrb,rhovan,drhovan,rhog,rhor,drhog,drhor)
       USE ions_base,                ONLY: na, nsp, nat
       USE cvan,                     ONLY: nvb
       USE uspp_param,               ONLY: nhm, nh
-      USE grid_dimensions,          ONLY: nr1, nr2, nr3, nr1x, nr2x, nr3x, nnr => nrxx
+      USE grid_dimensions,          ONLY: nr1, nr2, nr3, nr1x, nr2x, nr3x, nrxx
       USE electrons_base,           ONLY: nspin
       USE gvecb,                    ONLY: ngb, npb, nmb
       USE gvecp,                    ONLY: ngm
-      USE smallbox_grid_dimensions, ONLY: nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx, nnrb => nnrbx
+      USE smallbox_grid_dimensions, ONLY: nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx, nnrbx
       USE cell_base,                ONLY: ainv
       USE qgb_mod,                  ONLY: qgb
       USE dqgb_mod,                 ONLY: dqgb
@@ -783,12 +783,12 @@ SUBROUTINE drhov(irb,eigrb,rhovan,drhovan,rhog,rhor,drhog,drhor)
       IMPLICIT NONE
 ! input
       INTEGER,     INTENT(IN) ::  irb(3,nat)
-      REAL(DP),    INTENT(IN) ::  rhor(nnr,nspin)
+      REAL(DP),    INTENT(IN) ::  rhor(nrxx,nspin)
       REAL(DP),    INTENT(IN) ::  rhovan(nhm*(nhm+1)/2,nat,nspin)
       REAL(DP),    INTENT(IN) ::  drhovan(nhm*(nhm+1)/2,nat,nspin,3,3)
       COMPLEX(DP), INTENT(IN) ::  eigrb(ngb,nat), rhog(ngm,nspin)
 ! output
-      REAL(DP),    INTENT(OUT) :: drhor(nnr,nspin,3,3)
+      REAL(DP),    INTENT(OUT) :: drhor(nrxx,nspin,3,3)
       COMPLEX(DP), INTENT(OUT) :: drhog(ngm,nspin,3,3)
 ! local
       INTEGER i, j, isup, isdw, nfft, ifft, iv, jv, ig, ijv, is, iss,   &
@@ -808,7 +808,7 @@ SUBROUTINE drhov(irb,eigrb,rhovan,drhovan,rhog,rhor,drhog,drhor)
       DO j=1,3
          DO i=1,3
             DO iss=1,nspin
-               DO ir=1,nnr
+               DO ir=1,nrxx
                   drhor(ir,iss,i,j)=-rhor(ir,iss)*ainv(j,i)
                END DO
                DO ig=1,ngm
@@ -820,7 +820,7 @@ SUBROUTINE drhov(irb,eigrb,rhovan,drhovan,rhog,rhor,drhog,drhor)
 
       IF ( nvb < 0 ) RETURN
 
-      ALLOCATE( v( nnr ) )
+      ALLOCATE( v( nrxx ) )
 
       ci =( 0.0d0, 1.0d0 )
 
@@ -834,12 +834,12 @@ SUBROUTINE drhov(irb,eigrb,rhovan,drhovan,rhog,rhor,drhog,drhor)
                v(:) = (0.d0, 0.d0)
 
 !$omp parallel default(none) &
-!$omp          shared(nvb, na, nnrb, ngb, nh, eigrb, dfftb, irb, v, &
+!$omp          shared(nvb, na, nnrbx, ngb, nh, eigrb, dfftb, irb, v, &
 !$omp                 nmb, ci, npb, i, j, dqgb, qgb, nhm, rhovan, drhovan ) &
 !$omp          private(mytid, ntids, is, ia, nfft, ifft, iv, jv, ijv, ig, iss, isa, &
 !$omp                  qv, itid, dqgbt, dsumt, asumt )
 
-               ALLOCATE( qv( nnrb ) )
+               ALLOCATE( qv( nnrbx ) )
                ALLOCATE( dqgbt( ngb, 2 ) )
 
 #ifdef __OPENMP
@@ -937,7 +937,7 @@ SUBROUTINE drhov(irb,eigrb,rhovan,drhovan,rhog,rhor,drhog,drhor)
 
                iss = 1
 
-               DO ir=1,nnr
+               DO ir=1,nrxx
                   drhor(ir,iss,i,j) = drhor(ir,iss,i,j) + DBLE(v(ir))
                END DO
 !
@@ -959,7 +959,7 @@ SUBROUTINE drhov(irb,eigrb,rhovan,drhovan,rhog,rhor,drhog,drhor)
          DO i=1,3
             DO j=1,3
                v(:) = (0.d0, 0.d0)
-               ALLOCATE( qv( nnrb ) )
+               ALLOCATE( qv( nnrbx ) )
                ALLOCATE( dqgbt( ngb, 2 ) )
                isa=1
                DO is=1,nvb
@@ -1014,7 +1014,7 @@ SUBROUTINE drhov(irb,eigrb,rhovan,drhovan,rhog,rhor,drhog,drhor)
                DEALLOCATE( dqgbt )
                DEALLOCATE( qv )
 !
-               DO ir=1,nnr
+               DO ir=1,nrxx
                   drhor(ir,isup,i,j) = drhor(ir,isup,i,j) + DBLE(v(ir))
                   drhor(ir,isdw,i,j) = drhor(ir,isdw,i,j) +AIMAG(v(ir))
                ENDDO
@@ -1059,13 +1059,13 @@ SUBROUTINE rhov(irb,eigrb,rhovan,rhog,rhor)
       USE cvan,                     ONLY: nvb
       USE uspp_param,               ONLY: nh, nhm
       USE uspp,                     ONLY: deeq
-      USE grid_dimensions,          ONLY: nr1, nr2, nr3, nr1x, nr2x, nr3x, nnr => nrxx
+      USE grid_dimensions,          ONLY: nr1, nr2, nr3, nr1x, nr2x, nr3x, nrxx
       USE electrons_base,           ONLY: nspin
       USE gvecb,                    ONLY: npb, nmb, ngb
       USE gvecp,                    ONLY: ngm
       USE cell_base,                ONLY: omega, ainv
       USE small_box,                ONLY: omegab
-      USE smallbox_grid_dimensions, ONLY: nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx, nnrb => nnrbx
+      USE smallbox_grid_dimensions, ONLY: nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx, nnrbx
       USE control_flags,            ONLY: iprint, iprsta, tpre
       USE qgb_mod,                  ONLY: qgb
       USE recvecs_indexes,          ONLY: np, nm
@@ -1078,7 +1078,7 @@ SUBROUTINE rhov(irb,eigrb,rhovan,rhog,rhor)
       INTEGER,     INTENT(in) :: irb(3,nat)
       COMPLEX(DP), INTENT(in):: eigrb(ngb,nat)
       ! 
-      REAL(DP),     INTENT(inout):: rhor(nnr,nspin)
+      REAL(DP),     INTENT(inout):: rhor(nrxx,nspin)
       COMPLEX(DP),  INTENT(inout):: rhog(ngm,nspin)
 !
       INTEGER     :: isup, isdw, nfft, ifft, iv, jv, ig, ijv, is, iss, isa, ia, ir, i, j
@@ -1102,7 +1102,7 @@ SUBROUTINE rhov(irb,eigrb,rhovan,rhog,rhor)
       ci=(0.d0,1.d0)
 !
 !
-      ALLOCATE( v( nnr ) )
+      ALLOCATE( v( nrxx ) )
 
       ! private variable need to be initialized, otherwise
       ! outside the parallel region they have an undetermined value
@@ -1121,10 +1121,10 @@ SUBROUTINE rhov(irb,eigrb,rhovan,rhog,rhor)
          !
 
 !$omp parallel default(none) &
-!$omp          shared(nvb, na, nnrb, ngb, nh, rhovan, qgb, eigrb, dfftb, iprsta, omegab, irb, v, nr1b, &
+!$omp          shared(nvb, na, nnrbx, ngb, nh, rhovan, qgb, eigrb, dfftb, iprsta, omegab, irb, v, nr1b, &
 !$omp                 nr2b, nr3b, nmb, stdout, ci, npb, rhor ) &
 !$omp          private(mytid, ntids, is, ia, nfft, ifft, iv, jv, ijv, sumrho, qgbt, ig, iss, isa, ca, &
-!$omp                  qv, itid, ir, nnr )
+!$omp                  qv, itid, ir, nrxx )
 
          iss=1
          isa=1
@@ -1140,7 +1140,7 @@ SUBROUTINE rhov(irb,eigrb,rhovan,rhog,rhor)
 #endif
 
          ALLOCATE( qgbt( ngb, 2 ) )
-         ALLOCATE( qv( nnrb ) )
+         ALLOCATE( qv( nnrbx ) )
 
 
          DO is = 1, nvb
@@ -1243,7 +1243,7 @@ SUBROUTINE rhov(irb,eigrb,rhovan,rhog,rhor)
 
          iss = 1
 
-         DO ir=1,nnr
+         DO ir=1,nrxx
             rhor(ir,iss)=rhor(ir,iss)+DBLE(v(ir))        
          END DO
 
@@ -1285,7 +1285,7 @@ SUBROUTINE rhov(irb,eigrb,rhovan,rhog,rhor)
          v (:) = (0.d0, 0.d0)
 
          ALLOCATE( qgbt( ngb, 2 ) )
-         ALLOCATE( qv( nnrb ) )
+         ALLOCATE( qv( nnrbx ) )
 
          isa=1
          DO is=1,nvb
@@ -1342,7 +1342,7 @@ SUBROUTINE rhov(irb,eigrb,rhovan,rhog,rhor)
             END DO
          END DO
 !
-         DO ir=1,nnr
+         DO ir=1,nrxx
             rhor(ir,isup)=rhor(ir,isup)+DBLE(v(ir)) 
             rhor(ir,isdw)=rhor(ir,isdw)+AIMAG(v(ir)) 
          END DO
diff --git a/CPV/cplib.f90 b/CPV/cplib.f90
index e06abaeae..f4095bfb4 100644
--- a/CPV/cplib.f90
+++ b/CPV/cplib.f90
@@ -881,11 +881,11 @@ END FUNCTION
       USE cvan,                     ONLY: nvb
       USE ions_base,                ONLY: nat, nsp, na
       USE constants,                ONLY: pi, fpi
-      USE grid_dimensions,          ONLY: nr3, nnr => nrxx
+      USE grid_dimensions,          ONLY: nr3, nrxx
       USE gvecb,                    ONLY: ngb, npb, nmb, gxb
       USE small_box,                ONLY: omegab, tpibab
       USE smallbox_grid_dimensions, ONLY: nr1b, nr2b, nr3b, &
-                                          nr1bx, nr2bx, nr3bx, nnrb => nnrbx
+                                          nr1bx, nr2bx, nr3bx, nnrbx
       USE qgb_mod,                  ONLY: qgb
       USE electrons_base,           ONLY: nspin
       USE control_flags,            ONLY: iprint, thdyn, tfor, tprnfor
@@ -899,7 +899,7 @@ END FUNCTION
       INTEGER irb(3,nat)
       REAL(DP) rhovan(nhm*(nhm+1)/2,nat,nspin)
       COMPLEX(DP) eigrb(ngb,nat)
-      REAL(DP)  vr(nnr,nspin)
+      REAL(DP)  vr(nrxx,nspin)
 ! output
       REAL(DP)  fion(3,nat)
 ! local
@@ -924,7 +924,7 @@ END FUNCTION
 
 
 !$omp parallel default(none) &
-!$omp          shared(nvb, na, nnrb, ngb, nh, qgb, eigrb, dfftb, irb, vr, nmb, npb, ci, deeq, &
+!$omp          shared(nvb, na, nnrbx, ngb, nh, qgb, eigrb, dfftb, irb, vr, nmb, npb, ci, deeq, &
 !$omp                 fac, nspin ) &
 !$omp          private(mytid, ntids, is, ia, nfft, iv, jv, ijv, ig, isa, qv, itid, res, iss )
 
@@ -937,7 +937,7 @@ END FUNCTION
 #endif
 
 
-      ALLOCATE( qv( nnrb ) )
+      ALLOCATE( qv( nnrbx ) )
 !
 ! calculation of deeq_i,lm = \int V_eff(r) q_i,lm(r) dr
 !
@@ -1029,13 +1029,13 @@ END FUNCTION
          !     -----------------------------------------------------------------
 
 !$omp parallel default(none) &
-!$omp          shared(nvb, na, nnrb, ngb, nh, qgb, eigrb, dfftb, irb, vr, nmb, npb, ci, deeq, &
+!$omp          shared(nvb, na, nnrbx, ngb, nh, qgb, eigrb, dfftb, irb, vr, nmb, npb, ci, deeq, &
 !$omp                 fac, nspin, rhovan, tpibab, gxb, fvan ) &
 !$omp          private(mytid, ntids, is, ia, ik, nfft, iv, jv, ijv, ig, isa, qv, itid, res, iss, &
 !$omp                  fac1, fac2, facg1, facg2 )
 
 
-         ALLOCATE( qv( nnrb ) )
+         ALLOCATE( qv( nnrbx ) )
 
          iss=1
          isa=1
@@ -1134,7 +1134,7 @@ END FUNCTION
          !     case nspin=2: up and down spin fft's combined into a single fft
          !     -----------------------------------------------------------------
 
-         ALLOCATE( qv( nnrb ) )
+         ALLOCATE( qv( nnrbx ) )
          isup=1
          isdw=2
          isa=1
@@ -1527,8 +1527,7 @@ END FUNCTION
       USE mp, ONLY: mp_sum
       USE gvecw, ONLY: ngw
       USE reciprocal_vectors, ONLY: gstart
-      USE grid_dimensions, ONLY: nr1, nr2, nr3, &
-            nnr => nrxx
+      USE grid_dimensions, ONLY: nr1, nr2, nr3, nrxx
       USE cell_base, ONLY: omega
       USE cvan, ONLY: nvb, ish
       USE uspp, ONLY: nhsa => nkb, nhsavb=>nkbus, qq
@@ -1537,7 +1536,7 @@ END FUNCTION
 !
       IMPLICIT NONE
 ! input
-      REAL(dp) bec(nhsa,n), rhor(nnr,nspin)
+      REAL(dp) bec(nhsa,n), rhor(nrxx,nspin)
       COMPLEX(dp) c(ngw,nx)
 ! local variables
       INTEGER nup, ndw, ir, i, j, jj, ig, ia, is, iv, jv, inl, jnl
@@ -1576,7 +1575,7 @@ END FUNCTION
 !     Becke's formula for spin polarization
 !
       spin1 = 0.0d0
-      DO ir=1,nnr
+      DO ir=1,nrxx
          spin1 = spin1 - MIN(rhor(ir,1),rhor(ir,2))
       END DO
       CALL mp_sum( spin1, intra_image_comm )
@@ -1671,7 +1670,7 @@ END FUNCTION
       USE cell_base,          ONLY: a1, a2, a3, tpiba2, h, ainv
       USE reciprocal_vectors, ONLY: gstart, gg, g
       USE recvecs_indexes,    ONLY: np, nm
-      USE grid_dimensions,    ONLY: nr1, nr2, nr3, nnr => nrxx
+      USE grid_dimensions,    ONLY: nr1, nr2, nr3, nrxx
       USE smooth_grid_dimensions, ONLY: nrxxs
       USE electrons_base,   ONLY: nspin
       USE constants,        ONLY: pi, fpi, au_gpa
@@ -1700,8 +1699,8 @@ END FUNCTION
 !
       LOGICAL :: tlast, tfirst
       INTEGER :: nfi
-      REAL(DP)  rhor(nnr,nspin), rhos(nrxxs,nspin), fion(3,nat)
-      REAL(DP)  rhoc(nnr), tau0(3,nat)
+      REAL(DP)  rhor(nrxx,nspin), rhos(nrxxs,nspin), fion(3,nat)
+      REAL(DP)  rhoc(nrxx), tau0(3,nat)
       COMPLEX(DP) ei1(-nr1:nr1,nat), ei2(-nr2:nr2,nat),     &
      &                ei3(-nr3:nr3,nat), eigrb(ngb,nat),        &
      &                rhog(ngm,nspin), sfac(ngms,nsp)
@@ -1749,7 +1748,7 @@ END FUNCTION
       !
       ht = TRANSPOSE( h )
       !
-      ALLOCATE( v( nnr ) )
+      ALLOCATE( v( nrxx ) )
       ALLOCATE( vs( nrxxs ) )
       ALLOCATE( vtemp( ngm ) )
       ALLOCATE( rhotmp( ngm ) )
@@ -1960,12 +1959,12 @@ END FUNCTION
          iss = 1
          if (abivol.or.abisur) then
 !$omp parallel do
-            do ir=1,nnr
+            do ir=1,nrxx
                v(ir)=CMPLX( rhor( ir, iss ) + v_vol( ir ), 0.d0 ,kind=DP)
             end do           
          else
 !$omp parallel do
-            do ir=1,nnr
+            do ir=1,nrxx
                v(ir)=CMPLX( rhor( ir, iss ), 0.d0 ,kind=DP)
             end do
          end if
@@ -1987,13 +1986,13 @@ END FUNCTION
          isdw=2
          if (abivol.or.abisur) then
 !$omp parallel do
-            do ir=1,nnr
+            do ir=1,nrxx
                v(ir)=CMPLX ( rhor(ir,isup)+v_vol(ir), &
                              rhor(ir,isdw)+v_vol(ir),kind=DP)
             end do
          else
 !$omp parallel do
-            do ir=1,nnr
+            do ir=1,nrxx
                v(ir)=CMPLX (rhor(ir,isup),rhor(ir,isdw),kind=DP)
             end do
          end if
@@ -2050,7 +2049,7 @@ END FUNCTION
          CALL invfft('Dense',v, dfftp )
 !
 !$omp parallel do
-         DO ir=1,nnr
+         DO ir=1,nrxx
             rhor(ir,iss)=DBLE(v(ir))
          END DO
 !
@@ -2068,7 +2067,7 @@ END FUNCTION
 !
          CALL invfft('Dense',v, dfftp )
 !$omp parallel do
-         DO ir=1,nnr
+         DO ir=1,nrxx
             rhor(ir,isup)= DBLE(v(ir))
             rhor(ir,isdw)=AIMAG(v(ir))
          END DO
diff --git a/CPV/cplib_meta.f90 b/CPV/cplib_meta.f90
index 91a4bbefc..4fa9eff65 100644
--- a/CPV/cplib_meta.f90
+++ b/CPV/cplib_meta.f90
@@ -16,7 +16,7 @@
       use reciprocal_vectors,     only : g
       use gvecs,                  only : ngms, nms, nps
       use gvecw,                  only : ngw
-      use smooth_grid_dimensions, only : nnrs => nrxxs
+      use smooth_grid_dimensions, only : nrxxs
       use cell_base,              only : tpiba2
       USE metagga,                ONLY : kedtaus
       USE fft_interfaces,         ONLY : fwfft, invfft
@@ -24,7 +24,7 @@
 !
       implicit none
 !
-      complex(dp) c(ngw), ca(ngw), df(ngw), da(ngw),psi(nnrs)
+      complex(dp) c(ngw), ca(ngw), df(ngw), da(ngw),psi(nrxxs)
       integer iss1, iss2
       real(dp) fi, fip
 ! local variables
@@ -42,7 +42,7 @@
             end do
             call invfft('Wave',psi,dffts )
 !           on smooth grids--> grids for charge density
-            do ir=1, nnrs
+            do ir=1, nrxxs
                psi(ir) = CMPLX (kedtaus(ir,iss1)*DBLE(psi(ir)), &
                                 kedtaus(ir,iss2)*AIMAG(psi(ir)),kind=DP)
             end do
@@ -72,8 +72,7 @@
       use gvecw, only: ngw
       use reciprocal_vectors, only: g
       use recvecs_indexes, only: np, nm
-      use grid_dimensions, only: nr1, nr2, nr3, &
-            nr1x, nr2x, nr3x, nnr => nrxx
+      use grid_dimensions, only: nr1, nr2, nr3, nr1x, nr2x, nr3x, nrxx
       use cell_base
       use smooth_grid_dimensions, only: nrxxs
       use electrons_base, only: nx => nbspx, n => nbsp, f, ispin, nspin
@@ -213,7 +212,7 @@
          psi(nm(1:ngms))=CONJG(kedtaug(1:ngms,iss))
          psi(np(1:ngms))=      kedtaug(1:ngms,iss)
          call invfft('Dense',psi, dfftp )
-         kedtaur(1:nnr,iss)=DBLE(psi(1:nnr))
+         kedtaur(1:nrxx,iss)=DBLE(psi(1:nrxx))
 
       else 
 !     ==================================================================
@@ -229,8 +228,8 @@
             psi(np(ig))=kedtaug(ig,isup)+ci*kedtaug(ig,isdw)
          end do
          call invfft('Dense',psi, dfftp )
-         kedtaur(1:nnr,isup)= DBLE(psi(1:nnr))
-         kedtaur(1:nnr,isdw)=AIMAG(psi(1:nnr))
+         kedtaur(1:nrxx,isup)= DBLE(psi(1:nrxx))
+         kedtaur(1:nrxx,isdw)=AIMAG(psi(1:nrxx))
 
       endif
 
@@ -262,9 +261,8 @@
       use cell_base, only: omega
       use cell_base, only: a1, a2, a3, tpiba2
       use recvecs_indexes, only: np, nm
-      use grid_dimensions, only: nr1, nr2, nr3, &
-            nr1x, nr2x, nr3x, nnr => nrxx
-      use smooth_grid_dimensions, only: nr1s, nr2s, nr3s, nnrs => nrxxs
+      use grid_dimensions, only: nr1, nr2, nr3, nr1x, nr2x, nr3x, nrxx
+      use smooth_grid_dimensions, only: nr1s, nr2s, nr3s, nrxxs
       use electrons_base, only: nspin
       use constants, only: pi, fpi
       use energies, only: etot, eself, enl, ekin, epseu, esr, eht, exc
@@ -284,7 +282,7 @@
       integer iss, isup, isdw, ig, ir,i,j,k,is, ia
       real(dp) dkedxc(3,3) !metagga
       complex(dp)  fp, fm, ci
-      complex(dp)  v(nnr), vs(nnrs)
+      complex(dp)  v(nrxx), vs(nrxxs)
 !
       ci=(0.d0,1.d0)
 
@@ -309,7 +307,7 @@
 !
       if(nspin.eq.1) then
          iss=1
-         do ir=1,nnr
+         do ir=1,nrxx
             v(ir)=CMPLX(kedtaur(ir,iss),0.0d0,kind=DP)
          end do
          call fwfft('Dense',v, dfftp )
@@ -321,7 +319,7 @@
          isup=1
          isdw=2
 
-         v(1:nnr)=CMPLX(kedtaur(1:nnr,isup),kedtaur(1:nnr,isdw),kind=DP)
+         v(1:nrxx)=CMPLX(kedtaur(1:nrxx,isup),kedtaur(1:nrxx,isdw),kind=DP)
          call fwfft('Dense',v, dfftp )
          do ig=1,ngm
             fp=v(np(ig))+v(nm(ig))
@@ -342,7 +340,7 @@
 !
          call invfft('Smooth',vs, dffts )
 !
-         kedtaus(1:nnrs,iss)=DBLE(vs(1:nnrs))
+         kedtaus(1:nrxxs,iss)=DBLE(vs(1:nrxxs))
       else
          isup=1
          isdw=2
@@ -351,8 +349,8 @@
             vs(nms(ig))=CONJG(kedtaug(ig,isup)) +ci*conjg(kedtaug(ig,isdw))
          end do
          call invfft('Smooth',vs, dffts )
-         kedtaus(1:nnrs,isup)= DBLE(vs(1:nnrs))
-         kedtaus(1:nnrs,isdw)=AIMAG(vs(1:nnrs))
+         kedtaus(1:nrxxs,isup)= DBLE(vs(1:nrxxs))
+         kedtaus(1:nrxxs,isdw)=AIMAG(vs(1:nrxxs))
       endif
       !calculate dkedxc in real space on smooth grids  !metagga
       if(tpre) then
@@ -360,7 +358,7 @@
             do j=1,3
                do i=1,3
                   dkedxc(i,j)=0.d0
-                  do ir=1,nnrs
+                  do ir=1,nrxxs
                      !2.d0 : because kedtau = 0.5d0 d_Exc/d_kedtau
                       dkedxc(i,j)= dkedxc(i,j)+kedtaus(ir,iss)*2.d0*&
                            dkedtaus(ir,i,j,iss)
diff --git a/CPV/exch_corr.f90 b/CPV/exch_corr.f90
index 16388bba3..83b0a6a66 100644
--- a/CPV/exch_corr.f90
+++ b/CPV/exch_corr.f90
@@ -17,7 +17,7 @@
       use funct,           only : dft_is_gradient, dft_is_meta
       use gvecp,           only : ngm
       use gvecs,           only : ngms
-      use grid_dimensions, only : nr1, nr2, nr3, nnr => nrxx
+      use grid_dimensions, only : nr1, nr2, nr3, nrxx
       use cell_base,       only : ainv, omega, h
       use ions_base,       only : nsp
       use control_flags,   only : tpre, iprsta
@@ -47,7 +47,7 @@
       ! output
       ! rhor contains the exchange-correlation potential
       !
-      real(DP) :: rhor( nnr, nspin ), rhoc( nnr )
+      real(DP) :: rhor( nrxx, nspin ), rhoc( nrxx )
       real(DP) :: dxc( 3, 3 ), exc
       real(DP) :: dcc( 3, 3 ), drc( 3, 3 )
       !
@@ -68,7 +68,7 @@
       !
       if ( dft_is_gradient() ) then
          !
-         allocate( gradr( nnr, 3, nspin ) )
+         allocate( gradr( nrxx, 3, nspin ) )
          call fillgrad( nspin, rhog, gradr )
          ! 
       else
@@ -88,9 +88,9 @@
 
          !  allocate the sic_arrays
          !
-         ALLOCATE( self_rho( nnr, nspin ) )
+         ALLOCATE( self_rho( nrxx, nspin ) )
          ALLOCATE( self_rhog(ngm, nspin ) )
-         IF( dft_is_gradient() ) ALLOCATE( self_gradr( nnr, 3, nspin ) )
+         IF( dft_is_gradient() ) ALLOCATE( self_gradr( nrxx, 3, nspin ) )
 
          self_rho(:, 1) = rhor( :, 2)
          self_rho(:, 2) = rhor( :, 2)
@@ -109,14 +109,14 @@
 !
       if( dft_is_meta() ) then
          !
-         call tpssmeta( nnr, nspin, gradr, rhor, kedtaur, exc )
+         call tpssmeta( nrxx, nspin, gradr, rhor, kedtaur, exc )
          !
       else
          !
-         CALL exch_corr_cp(nnr, nspin, gradr, rhor, exc)
+         CALL exch_corr_cp(nrxx, nspin, gradr, rhor, exc)
          !
          IF ( ttsic ) THEN
-            CALL exch_corr_cp(nnr, nspin, self_gradr, self_rho, self_exc)
+            CALL exch_corr_cp(nrxx, nspin, self_gradr, self_rho, self_exc)
             self_exc = sic_alpha * (exc - self_exc)
             exc = exc - self_exc
          END IF
@@ -142,7 +142,7 @@
          do iss = 1, nspin
             do j=1,3
                do i=1,3
-                  do ir=1,nnr
+                  do ir=1,nrxx
                      dxc(i,j) = dxc(i,j) + rhor( ir, iss ) * drhor( ir, iss, i, j )
                   end do
                end do
@@ -220,7 +220,7 @@
          !
          dcc = 0.0d0
          !
-         IF( nlcc_any ) CALL  denlcc( nnr, nspin, rhor, sfac, drhocg, dcc )
+         IF( nlcc_any ) CALL  denlcc( nrxx, nspin, rhor, sfac, drhocg, dcc )
          !
          ! DEBUG
          !
@@ -234,7 +234,7 @@
             IF( nlcc_any ) THEN
                do j=1,3
                   do i=1,3
-                     do ir=1,nnr
+                     do ir=1,nrxx
                         drc(i,j) = drc(i,j) + rhor( ir, iss ) * rhoc( ir ) * ainv(j,i)
                      end do
                   end do
@@ -270,7 +270,7 @@
       use reciprocal_vectors, only: g
       use recvecs_indexes, only: np, nm
       use gvecp, only: ngm
-      use grid_dimensions, only: nr1, nr2, nr3, nnr => nrxx, nr1x, nr2x, nr3x
+      use grid_dimensions, only: nr1, nr2, nr3, nrxx, nr1x, nr2x, nr3x
       use cell_base, only: ainv, tpiba, omega
       use cp_main_variables, only: drhog
       USE fft_interfaces, ONLY: fwfft, invfft
@@ -279,7 +279,7 @@
       implicit none  
 ! input                   
       integer nspin
-      real(DP)    :: gradr( nnr, 3, nspin ), rhor( nnr, nspin ), dexc( 3, 3 )
+      real(DP)    :: gradr( nrxx, 3, nspin ), rhor( nrxx, nspin ), dexc( 3, 3 )
       complex(DP) :: rhog( ngm, nspin )
 !
       complex(DP), allocatable:: v(:)
@@ -287,7 +287,7 @@
       complex(DP) ::  ci, fp, fm
       integer :: iss, ig, ir, i,j
 !
-      allocate(v(nnr))
+      allocate(v(nrxx))
       allocate(x(ngm))
       allocate(vtemp(ngm))
       !
@@ -299,7 +299,7 @@
 !     _________________________________________________________________
 !     second part xc-potential: 3 forward ffts
 !
-         do ir=1,nnr
+         do ir=1,nrxx
             v(ir)=CMPLX(gradr(ir,1,iss),0.d0,kind=DP)
          end do
          call fwfft('Dense',v, dfftp )
@@ -320,7 +320,7 @@
             end do
          endif
 !
-         do ir=1,nnr
+         do ir=1,nrxx
             v(ir)=CMPLX(gradr(ir,2,iss),gradr(ir,3,iss),kind=DP)
          end do
          call fwfft('Dense',v, dfftp )
@@ -355,7 +355,7 @@
 !     _________________________________________________________________
 !     second part xc-potential: 1 inverse fft
 !
-         do ig=1,nnr
+         do ig=1,nrxx
             v(ig)=(0.0d0,0.0d0)
          end do
          do ig=1,ngm
@@ -363,7 +363,7 @@
             v(nm(ig))=CONJG(x(ig))
          end do
          call invfft('Dense',v, dfftp )
-         do ir=1,nnr
+         do ir=1,nrxx
             rhor(ir,iss)=rhor(ir,iss)-DBLE(v(ir))
          end do
       end do
@@ -383,17 +383,17 @@
 !
 !=----------------------------------------------------------------------------=!
 
-subroutine exch_corr_wrapper(nnr, nspin, grhor, rhor, etxc, v, h)
+subroutine exch_corr_wrapper(nrxx, nspin, grhor, rhor, etxc, v, h)
   use kinds, only: DP
   use funct, only: dft_is_gradient, get_igcc, &
                    xc, xc_spin, gcxc, gcx_spin, gcc_spin, gcc_spin_more
   implicit none
-  integer, intent(in) :: nnr
+  integer, intent(in) :: nrxx
   integer, intent(in) :: nspin
-  real(DP), intent(in) :: grhor( nnr, 3, nspin )
-  real(DP) :: h( nnr, nspin, nspin )
-  real(DP), intent(in) :: rhor( nnr, nspin )
-  real(DP) :: v( nnr, nspin )
+  real(DP), intent(in) :: grhor( nrxx, 3, nspin )
+  real(DP) :: h( nrxx, nspin, nspin )
+  real(DP), intent(in) :: rhor( nrxx, nspin )
+  real(DP) :: v( nrxx, nspin )
   real(DP) :: etxc
   integer :: ir, is, k
   real(DP) :: rup, rdw, ex, ec, vx(2), vc(2)
@@ -419,7 +419,7 @@ subroutine exch_corr_wrapper(nnr, nspin, grhor, rhor, etxc, v, h)
      ! spin-unpolarized case
      !
 !$omp parallel do private( rhox, arhox, ex, ec, vx, vc ), reduction(+:etxc)
-     do ir = 1, nnr
+     do ir = 1, nrxx
         rhox = rhor (ir, nspin)
         arhox = abs (rhox)
         if (arhox.gt.1.d-30) then
@@ -439,7 +439,7 @@ subroutine exch_corr_wrapper(nnr, nspin, grhor, rhor, etxc, v, h)
      neg (1) = 0
      neg (2) = 0
      neg (3) = 0
-     do ir = 1, nnr
+     do ir = 1, nrxx
         rhox = rhor(ir,1) + rhor(ir,2)
         arhox = abs(rhox)
         if (arhox.gt.1.d-30) then
@@ -466,7 +466,7 @@ subroutine exch_corr_wrapper(nnr, nspin, grhor, rhor, etxc, v, h)
 
   if( debug_xc ) then
     open(unit=17,form='unformatted')
-    write(17) nnr, nspin
+    write(17) nrxx, nspin
     write(17) rhor
     write(17) grhor
     close(17)
@@ -483,7 +483,7 @@ subroutine exch_corr_wrapper(nnr, nspin, grhor, rhor, etxc, v, h)
        !
 !$omp parallel do &
 !$omp private( is, grho2, arho, segno, sx, sc, v1x, v2x, v1c, v2c  ), reduction(+:etxc)
-       do k = 1, nnr
+       do k = 1, nrxx
           !
           grho2 (1) = grhor(k, 1, 1)**2 + grhor(k, 2, 1)**2 + grhor(k, 3, 1)**2
           arho = abs (rhor (k, 1) )
@@ -512,7 +512,7 @@ subroutine exch_corr_wrapper(nnr, nspin, grhor, rhor, etxc, v, h)
        !
        !    spin-polarised case
        !
-       do k = 1, nnr
+       do k = 1, nrxx
           do is = 1, nspin
              grho2 (is) = grhor(k, 1, is)**2 + grhor(k, 2, is)**2 + grhor(k, 3, is)**2
           enddo
@@ -583,28 +583,28 @@ end subroutine exch_corr_wrapper
 !
 !=----------------------------------------------------------------------------=!
 
-subroutine exch_corr_cp(nnr,nspin,grhor,rhor,etxc)
+subroutine exch_corr_cp(nrxx,nspin,grhor,rhor,etxc)
   use kinds, only: DP
   use funct, only: dft_is_gradient
   implicit none
-  integer, intent(in) :: nnr
+  integer, intent(in) :: nrxx
   integer, intent(in) :: nspin
-  real(DP) :: grhor( nnr, 3, nspin )
-  real(DP) :: rhor( nnr, nspin )
+  real(DP) :: grhor( nrxx, 3, nspin )
+  real(DP) :: rhor( nrxx, nspin )
   real(DP) :: etxc
   integer :: k, ipol
   real(DP) ::  grup, grdw
   real(DP), allocatable :: v(:,:)
   real(DP), allocatable :: h(:,:,:)
   !
-  allocate( v( nnr, nspin ) )
+  allocate( v( nrxx, nspin ) )
   if( dft_is_gradient() ) then
-    allocate( h( nnr, nspin, nspin ) )
+    allocate( h( nrxx, nspin, nspin ) )
   else
     allocate( h( 1, 1, 1 ) )
   endif
   !
-  call exch_corr_wrapper(nnr,nspin,grhor,rhor,etxc,v,h)
+  call exch_corr_wrapper(nrxx,nspin,grhor,rhor,etxc,v,h)
 
   if( dft_is_gradient() ) then
      !
@@ -615,7 +615,7 @@ subroutine exch_corr_cp(nnr,nspin,grhor,rhor,etxc)
         !
         do ipol = 1, 3
 !$omp do
-           do k = 1, nnr
+           do k = 1, nrxx
               grhor (k, ipol, 1) = h (k, 1, 1) * grhor (k, ipol, 1)
            enddo
 !$omp end do
@@ -626,7 +626,7 @@ subroutine exch_corr_cp(nnr,nspin,grhor,rhor,etxc)
         !
         do ipol = 1, 3
 !$omp do
-           do k = 1, nnr
+           do k = 1, nrxx
              grup = grhor (k, ipol, 1)
              grdw = grhor (k, ipol, 2)
              grhor (k, ipol, 1) = h (k, 1, 1) * grup + h (k, 1, 2) * grdw
diff --git a/CPV/fft.f90 b/CPV/fft.f90
index acf5d1f1d..727941e08 100644
--- a/CPV/fft.f90
+++ b/CPV/fft.f90
@@ -254,17 +254,16 @@
 ! nfft=2  add imaginary part of qv(r) to real part of array vr(r)
 !
       USE kinds, ONLY: dp
-      USE grid_dimensions, ONLY: nr1, nr2, nr3, &
-            nr1x, nr2x, nnr => nrxx
+      USE grid_dimensions, ONLY: nr1, nr2, nr3, nr1x, nr2x, nrxx
       USE smallbox_grid_dimensions, ONLY: nr1b, nr2b, nr3b, &
-            nr1bx, nr2bx, nnrb => nnrbx
+            nr1bx, nr2bx, nnrbx
       USE fft_base, ONLY: dfftp
       USE mp_global, ONLY: me_image
 
       IMPLICIT NONE
       INTEGER, INTENT(in):: nfft, irb(3)
-      REAL(dp), INTENT(in):: qv(2,nnrb)
-      COMPLEX(dp), INTENT(inout):: vr(nnr)
+      REAL(dp), INTENT(in):: qv(2,nnrbx)
+      COMPLEX(dp), INTENT(inout):: vr(nrxx)
 !
       INTEGER ir1, ir2, ir3, ir, ibig1, ibig2, ibig3, ibig
       INTEGER me
@@ -312,18 +311,17 @@
 ! irb   : position of the box in the dense grid
 !
       USE kinds, ONLY: dp
-      USE grid_dimensions, ONLY: nr1, nr2, nr3, &
-            nr1x, nr2x, nnr => nrxx
+      USE grid_dimensions, ONLY: nr1, nr2, nr3, nr1x, nr2x, nrxx
       USE smallbox_grid_dimensions, ONLY: nr1b, nr2b, nr3b, &
-            nr1bx, nr2bx, nnrb => nnrbx
+            nr1bx, nr2bx, nnrbx
       USE fft_base, ONLY: dfftp
       USE mp_global, ONLY: me_image
       !
       IMPLICIT NONE
       !
       INTEGER, INTENT(in):: irb(3)
-      COMPLEX(dp), INTENT(in):: qv(nnrb)
-      COMPLEX(dp), INTENT(inout):: v(nnr)
+      COMPLEX(dp), INTENT(in):: qv(nnrbx)
+      COMPLEX(dp), INTENT(inout):: v(nrxx)
 !
       INTEGER ir1, ir2, ir3, ir, ibig1, ibig2, ibig3, ibig
       INTEGER me
@@ -370,15 +368,14 @@
 ! Parallel execution: remember to sum the contributions from other nodes
 !
       USE kinds, ONLY: dp
-      USE grid_dimensions, ONLY: nr1, nr2, nr3, &
-            nr1x, nr2x, nnr => nrxx
+      USE grid_dimensions, ONLY: nr1, nr2, nr3, nr1x, nr2x, nrxx
       USE smallbox_grid_dimensions, ONLY: nr1b, nr2b, nr3b, &
-            nr1bx, nr2bx, nnrb => nnrbx
+            nr1bx, nr2bx, nnrbx
       USE fft_base, ONLY: dfftp
       USE mp_global, ONLY: me_image
       IMPLICIT NONE
       INTEGER, INTENT(in):: nfft, irb(3)
-      REAL(dp), INTENT(in):: qv(2,nnrb), vr(nnr)
+      REAL(dp), INTENT(in):: qv(2,nnrbx), vr(nrxx)
 !
       INTEGER ir1, ir2, ir3, ir, ibig1, ibig2, ibig3, ibig
       INTEGER me
diff --git a/CPV/gradrho.f90 b/CPV/gradrho.f90
index 6625cb352..55486f3ad 100644
--- a/CPV/gradrho.f90
+++ b/CPV/gradrho.f90
@@ -21,29 +21,26 @@ SUBROUTINE gradrho(nspin,rhog,drho,d2rho,dxdyrho,dxdzrho,dydzrho)
       use recvecs_indexes
       USE fft_interfaces, ONLY: invfft
       USE fft_base,       ONLY: dfftp
-      use grid_dimensions, only : nr1, nr2, nr3, nr1x, nr2x, nr3x,      &
-     &                            nnr=> nrxx
-!      use grid_dimensions, only: nr1, nr2, nr3, &
-!            nr1x, nr2x, nr3x, nnr => nrxx
+      use grid_dimensions, only : nr1, nr2, nr3, nr1x, nr2x, nr3x, nrxx
 !
       implicit none
 ! input
       integer nspin
       complex(kind=8) rhog(ngm,nspin)
 ! output
-      real(kind=8)    drho(3,nnr), d2rho(3,nnr),     &
-     &                dxdyrho(nnr), dxdzrho(nnr),    &
-     &                dydzrho(nnr)
+      real(kind=8)    drho(3,nrxx), d2rho(3,nrxx),     &
+     &                dxdyrho(nrxx), dxdzrho(nrxx),    &
+     &                dydzrho(nrxx)
 ! local
       complex(kind=8), allocatable:: v(:), w(:)
       complex(kind=8) ci
       integer iss, ig, ir, j
 !
 !
-      allocate(v(nnr))
-      allocate(w(nnr))
+      allocate(v(nrxx))
+      allocate(w(nrxx))
       ci=(0.0d0,1.0d0)
-      do ir = 1,nnr
+      do ir = 1,nrxx
          do j = 1,3
             drho(j,ir) = 0.d0
             d2rho(j,ir) = 0.d0
@@ -54,7 +51,7 @@ SUBROUTINE gradrho(nspin,rhog,drho,d2rho,dxdyrho,dxdzrho,dydzrho)
       end do
       do iss=1,nspin
 
-         do ig=1,nnr
+         do ig=1,nrxx
             v(ig)=(0.0d0,0.0d0)
             w(ig)=(0.0d0,0.0d0)
          end do
@@ -66,12 +63,12 @@ SUBROUTINE gradrho(nspin,rhog,drho,d2rho,dxdyrho,dxdzrho,dydzrho)
          end do
          call invfft('Dense',v, dfftp )
          call invfft('Dense',w, dfftp )
-         do ir=1,nnr
+         do ir=1,nrxx
             drho(1,ir)=drho(1,ir)+real(v(ir))
             d2rho(1,ir)=d2rho(1,ir)+real(w(ir))
          end do
 !
-         do ig=1,nnr
+         do ig=1,nrxx
             v(ig)=(0.0d0,0.0d0)
             w(ig)=(0.0d0,0.0d0)
          end do
@@ -87,14 +84,14 @@ SUBROUTINE gradrho(nspin,rhog,drho,d2rho,dxdyrho,dxdzrho,dydzrho)
          end do
          call invfft('Dense',v, dfftp )
          call invfft('Dense',w, dfftp )
-         do ir=1,nnr
+         do ir=1,nrxx
             drho(2,ir)=drho(2,ir)+real(v(ir))
             drho(3,ir)=drho(3,ir)+aimag(v(ir))
             d2rho(2,ir)=d2rho(2,ir)+real(w(ir))
             d2rho(3,ir)=d2rho(3,ir)+aimag(w(ir))
          end do
 
-         do ig=1,nnr
+         do ig=1,nrxx
             v(ig)=(0.0d0,0.0d0)
          end do
          do ig=1,ngm
@@ -102,11 +99,11 @@ SUBROUTINE gradrho(nspin,rhog,drho,d2rho,dxdyrho,dxdzrho,dydzrho)
             v(nm(ig))=conjg(v(np(ig)))
          end do
          call invfft('Dense',v, dfftp )
-         do ir=1,nnr
+         do ir=1,nrxx
             dxdyrho(ir)=dxdyrho(ir)+real(v(ir))
          end do
 !
-         do ig=1,nnr
+         do ig=1,nrxx
             v(ig)=(0.0d0,0.0d0)
          end do
          do ig=1,ngm
@@ -117,7 +114,7 @@ SUBROUTINE gradrho(nspin,rhog,drho,d2rho,dxdyrho,dxdzrho,dydzrho)
      &                        ci*conjg(g(2,ig)*g(3,ig)*rhog(ig,iss)))
          end do
          call invfft('Dense',v, dfftp )
-         do ir=1,nnr
+         do ir=1,nrxx
             dxdzrho(ir)=dxdzrho(ir)+real(v(ir))
             dydzrho(ir)=dydzrho(ir)+aimag(v(ir))
          end do
diff --git a/CPV/init.f90 b/CPV/init.f90
index b2f81d572..79e6df1b2 100644
--- a/CPV/init.f90
+++ b/CPV/init.f90
@@ -22,12 +22,12 @@
       USE kinds,                    ONLY: dp
       use io_global,                only: stdout, ionode
       use control_flags,            only: gamma_only
-      use grid_dimensions,          only: nr1, nr2, nr3, nr1x, nr2x, nr3x, nnr => nrxx
+      use grid_dimensions,          only: nr1, nr2, nr3, nr1x, nr2x, nr3x
       use cell_base,                only: ainv, a1, a2, a3
       use cell_base,                only: omega, alat
       use small_box,                only: a1b, a2b, a3b, omegab, ainvb, tpibab, small_box_set
       use small_box,                only: alatb, b1b, b2b, b3b
-      use smallbox_grid_dimensions, only: nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx, nnrb => nnrbx
+      use smallbox_grid_dimensions, only: nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx
       use smooth_grid_dimensions,   only: nr1s, nr2s, nr3s, nr1sx, nr2sx, nr3sx, nrxxs
       USE grid_subroutines,         ONLY: realspace_grids_init, realspace_grids_para
       USE reciprocal_vectors,       ONLY: mill_g, g2_g, bi1, bi2, bi3
diff --git a/CPV/inner_loop_cold.f90 b/CPV/inner_loop_cold.f90
index b636cbbb0..24c6a745a 100644
--- a/CPV/inner_loop_cold.f90
+++ b/CPV/inner_loop_cold.f90
@@ -45,14 +45,12 @@
       USE cvan,           ONLY: nvb, ish
       USE ions_base,      ONLY: na, nat, pmass, nax, nsp, rcmax
       USE grid_dimensions, &
-                          ONLY: nnr => nrxx, nr1, nr2, nr3
+                          ONLY: nrxx, nr1, nr2, nr3
       USE cell_base,      ONLY: ainv, a1, a2, a3
       USE cell_base,      ONLY: omega, alat
       USE cell_base,      ONLY: h, hold, deth, wmass, tpiba2
       USE smooth_grid_dimensions, &
                           ONLY: nrxxs, nr1s, nr2s, nr3s
-      USE smallbox_grid_dimensions, &
-                          ONLY: nnrb => nnrbx, nr1b, nr2b, nr3b
       USE local_pseudo,   ONLY: vps, rhops
       USE io_global,      ONLY: io_global_start, stdout, ionode, &
                                 ionode_id
@@ -91,11 +89,11 @@
 
       INTEGER                :: irb( 3, nat )
       COMPLEX (kind=DP)           :: eigrb( ngb, nat )
-      REAL(kind=DP)               :: rhor( nnr, nspin )
-      REAL(kind=DP)               :: vpot( nnr, nspin )
+      REAL(kind=DP)               :: rhor( nrxx, nspin )
+      REAL(kind=DP)               :: vpot( nrxx, nspin )
       COMPLEX(kind=DP)            :: rhog( ngm, nspin )
       REAL(kind=DP)               :: rhos( nrxxs, nspin )
-      REAL(kind=DP)               :: rhoc( nnr )
+      REAL(kind=DP)               :: rhoc( nrxx )
       COMPLEX(kind=DP)            :: ei1( nr1:nr1, nat )
       COMPLEX(kind=DP)            :: ei2( nr2:nr2, nat )
       COMPLEX(kind=DP)            :: ei3( nr3:nr3, nat )
@@ -359,14 +357,12 @@
       USE cvan,           ONLY: nvb, ish
       USE ions_base,      ONLY: na, nat, pmass, nax, nsp, rcmax
       USE grid_dimensions, &
-                          ONLY: nnr => nrxx, nr1, nr2, nr3
+                          ONLY: nrxx, nr1, nr2, nr3
       USE cell_base,      ONLY: ainv, a1, a2, a3
       USE cell_base,      ONLY: omega, alat
       USE cell_base,      ONLY: h, hold, deth, wmass, tpiba2
       USE smooth_grid_dimensions, &
                           ONLY: nrxxs, nr1s, nr2s, nr3s
-      USE smallbox_grid_dimensions, &
-                          ONLY: nnrb => nnrbx, nr1b, nr2b, nr3b
       USE local_pseudo,   ONLY: vps, rhops
       USE io_global,      ONLY: io_global_start, stdout, ionode, &
                                 ionode_id
@@ -398,11 +394,11 @@
 
       INTEGER                :: irb( 3, nat )
       COMPLEX (kind=DP)           :: eigrb( ngb, nat )
-      REAL(kind=DP)               :: rhor( nnr, nspin )
-      REAL(kind=DP)               :: vpot( nnr, nspin )
+      REAL(kind=DP)               :: rhor( nrxx, nspin )
+      REAL(kind=DP)               :: vpot( nrxx, nspin )
       COMPLEX(kind=DP)            :: rhog( ngm, nspin )
       REAL(kind=DP)               :: rhos( nrxxs, nspin )
-      REAL(kind=DP)               :: rhoc( nnr )
+      REAL(kind=DP)               :: rhoc( nrxx )
       COMPLEX(kind=DP)            :: ei1( nr1:nr1, nat )
       COMPLEX(kind=DP)            :: ei2( nr2:nr2, nat )
       COMPLEX(kind=DP)            :: ei3( nr3:nr3, nat )
@@ -557,14 +553,12 @@
       USE cvan,           ONLY: nvb, ish
       USE ions_base,      ONLY: na, nat, pmass, nax, nsp, rcmax
       USE grid_dimensions, &
-                          ONLY: nnr => nrxx, nr1, nr2, nr3
+                          ONLY: nrxx, nr1, nr2, nr3
       USE cell_base,      ONLY: ainv, a1, a2, a3
       USE cell_base,      ONLY: omega, alat
       USE cell_base,      ONLY: h, hold, deth, wmass, tpiba2
       USE smooth_grid_dimensions, &
                           ONLY: nrxxs, nr1s, nr2s, nr3s
-      USE smallbox_grid_dimensions, &
-                          ONLY: nnrb => nnrbx, nr1b, nr2b, nr3b
       USE local_pseudo,   ONLY: vps, rhops
       USE io_global,      ONLY: io_global_start, stdout, ionode, &
                                 ionode_id
diff --git a/CPV/nlcc.f90 b/CPV/nlcc.f90
index 1649df83d..f69056bff 100644
--- a/CPV/nlcc.f90
+++ b/CPV/nlcc.f90
@@ -185,7 +185,7 @@
       USE kinds,           ONLY: DP
       use electrons_base,  only: nspin
       use gvecb,           only: gxb, ngb, npb, nmb
-      use grid_dimensions, only: nr1, nr2, nr3, nnr => nrxx
+      use grid_dimensions, only: nr1, nr2, nr3, nrxx
       use cell_base,       only: omega
       use ions_base,       only: nsp, na, nat
       use small_box,       only: tpibab
@@ -195,14 +195,14 @@
       use fft_base,        only: dfftb
       use reciprocal_vectors, only: gstart
       use smallbox_grid_dimensions, only: nr1b, nr2b, nr3b, &
-            nr1bx, nr2bx, nr3bx, nnrb => nnrbx
+            nr1bx, nr2bx, nr3bx, nnrbx
 
       implicit none
 
 ! input
       integer, intent(in)        :: irb(3,nat)
       complex(dp), intent(in):: eigrb(ngb,nat)
-      real(dp), intent(in)   :: vxc(nnr,nspin)
+      real(dp), intent(in)   :: vxc(nrxx,nspin)
 ! output
       real(dp), intent(inout):: fion1(3,nat)
 ! local
@@ -224,13 +224,13 @@
       fac = omega/DBLE(nr1*nr2*nr3*nspin)
 
 !$omp parallel default(none) &      
-!$omp          shared(nsp, na, nnrb, ngb, eigrb, dfftb, irb, nmb, npb, ci, rhocb, &
+!$omp          shared(nsp, na, nnrbx,ngb, eigrb, dfftb, irb, nmb, npb, ci, rhocb, &
 !$omp                 gxb, nat, fac, upf, vxc, nspin, tpibab, fion1 ) &
 !$omp          private(mytid, ntids, is, ia, nfft, ig, isa, qv, itid, res, ix, fcc, facg, iss )
 
 
       allocate( fcc( 3, nat ) )
-      allocate( qv( nnrb ) )
+      allocate( qv( nnrbx ) )
 
       fcc(:,:) = 0.d0
 
@@ -338,21 +338,21 @@
       use kinds, only: dp
       use ions_base,       only: nsp, na, nat
       use uspp_param,      only: upf
-      use grid_dimensions, only: nr3, nnr => nrxx
+      use grid_dimensions, only: nr3, nrxx
       use gvecb,           only: ngb, npb, nmb
       use control_flags,   only: iprint
       use core,            only: rhocb
       use fft_interfaces,  only: invfft
       use fft_base,        only: dfftb
       use smallbox_grid_dimensions, only: nr1b, nr2b, nr3b, &
-            nr1bx, nr2bx, nr3bx, nnrb => nnrbx
+            nr1bx, nr2bx, nr3bx, nnrbx
 
       implicit none
 ! input
       integer, intent(in)        :: irb(3,nat)
       complex(dp), intent(in):: eigrb(ngb,nat)
 ! output
-      real(dp), intent(out)  :: rhoc(nnr)
+      real(dp), intent(out)  :: rhoc(nrxx)
 ! local
       integer nfft, ig, is, ia, isa
       complex(dp) ci
@@ -367,15 +367,15 @@
       call start_clock( 'set_cc' )
       ci=(0.d0,1.d0)
 
-      allocate( wrk1 ( nnr ) )
+      allocate( wrk1 ( nrxx ) )
       wrk1 (:) = (0.d0, 0.d0)
 !
 !$omp parallel default(none) &      
-!$omp          shared(nsp, na, nnrb, ngb, eigrb, dfftb, irb, nmb, npb, ci, rhocb, &
+!$omp          shared(nsp, na, nnrbx, ngb, eigrb, dfftb, irb, nmb, npb, ci, rhocb, &
 !$omp                 nat, upf, wrk1 ) &
 !$omp          private(mytid, ntids, is, ia, nfft, ig, isa, qv, itid )
 
-      allocate( qv ( nnrb ) )
+      allocate( qv ( nnrbx ) )
 !
       isa = 0
 
@@ -442,7 +442,7 @@
 
 !$omp end parallel
 
-      call dcopy( nnr, wrk1, 2, rhoc, 1 )
+      call dcopy( nrxx, wrk1, 2, rhoc, 1 )
 
       deallocate( wrk1 )
 !
diff --git a/CPV/vanderwaals.f90 b/CPV/vanderwaals.f90
index b4d777ccc..86db3a3ad 100644
--- a/CPV/vanderwaals.f90
+++ b/CPV/vanderwaals.f90
@@ -23,7 +23,7 @@
 !---------------------------------
       subroutine VdW(evdw, taus, nat, na, nsp, fion, box) 
 
-      USE constants, ONLY: au => BOHR_RADIUS_ANGS
+      USE constants, ONLY: bohr_radius_angs
       USE cell_base, ONLY: s_to_r, boxdimensions
       USE mp_global, ONLY: me_image, root_image
 
@@ -130,13 +130,13 @@
 
                     ff = alp * (rcc - dist) 
                     dist6 = dist**6
-                    fun = - CSP1 / dist6 * cutofun_vdw(ff) / (au)**6
+                    fun = - CSP1 / dist6 * cutofun_vdw(ff) / bohr_radius_angs**6
 
                     if(dist.lt.rcut) then
                       fact = 0.d0 
                     else
                       fact = (6.d0 * CSP1/dist**7 * cutofun_vdw(ff) + &
-                      alp * dcutofun_vdw(ff) * CSP1/dist6) / (au)**6
+                      alp * dcutofun_vdw(ff) * CSP1/dist6) / bohr_radius_angs**6
                     endif
 
                     evdw = evdw + fun
diff --git a/CPV/vol_clu.f90 b/CPV/vol_clu.f90
index 96eec0d72..14f0edd2c 100644
--- a/CPV/vol_clu.f90
+++ b/CPV/vol_clu.f90
@@ -29,7 +29,7 @@ SUBROUTINE vol_clu(rho_real,rho_g,s_fac,flag)
       use control_flags,  only: tpre
       use fft_base,       ONLY : dfftp
       USE fft_interfaces, ONLY: invfft
-      use grid_dimensions,only: nr1, nr2, nr3, nr1x, nr2x, nr3x, nnr => nrxx
+      use grid_dimensions,only: nr1, nr2, nr3, nr1x, nr2x, nr3x, nrxx
       use pres_ai_mod, only: rho_thr, n_cntr, cntr, step_rad, fill_vac, &
      &                       delta_eps, delta_sigma, axis,              &
      &                       abisur, dthr, Surf_t, rho_gaus, v_vol,     &
@@ -52,7 +52,7 @@ SUBROUTINE vol_clu(rho_real,rho_g,s_fac,flag)
       real(kind=8) dx, dxx, xcc(4500)
       real(kind=8) weight0, wpiu, wmeno, maxr, minr
       real(kind=8) tau00(3), dist
-      real(kind=8) rho_real(nnr,nspin), rhoc
+      real(kind=8) rho_real(nrxx,nspin), rhoc
       real(kind=8) alfa(nsx), alfa0, sigma, hgt 
       real(kind=8) pos_cry(3), pos_car(3), pos_aux(3)
       real(kind=8) pos_cry0(3), dpvdh(3,3)
@@ -79,12 +79,12 @@ SUBROUTINE vol_clu(rho_real,rho_g,s_fac,flag)
       integer shift(nproc), incr(nproc),  ppp(nproc) 
       integer displs(nproc), ip, me
 #endif
-      if (abisur) allocate(drho(3,nnr))
-      if (abisur) allocate(d2rho(3,nnr))
-      if (abisur) allocate(dxdyrho(nnr))
-      if (abisur) allocate(dxdzrho(nnr))
-      if (abisur) allocate(dydzrho(nnr))
-      allocate(psi(nnr))
+      if (abisur) allocate(drho(3,nrxx))
+      if (abisur) allocate(d2rho(3,nrxx))
+      if (abisur) allocate(dxdyrho(nrxx))
+      if (abisur) allocate(dxdzrho(nrxx))
+      if (abisur) allocate(dydzrho(nrxx))
+      allocate(psi(nrxx))
 
       call start_clock( 'vol_clu' )
 
@@ -259,7 +259,7 @@ SUBROUTINE vol_clu(rho_real,rho_g,s_fac,flag)
             psi(nm(ig)) = conjg(rhotmp(ig,1))
          end do
          call invfft('Dense',psi, dfftp )
-         do ir = 1,nnr
+         do ir = 1,nrxx
             rho_gaus(ir) = real(psi(ir))
          end do
       else            
@@ -268,7 +268,7 @@ SUBROUTINE vol_clu(rho_real,rho_g,s_fac,flag)
             psi(nm(ig)) = conjg(rhotmp(ig,1)) + ci*conjg(rhotmp(ig,2))
          end do
          call invfft('Dense',psi, dfftp )
-         do ir = 1,nnr
+         do ir = 1,nrxx
             rho_gaus(ir) = real(psi(ir))+aimag(psi(ir))
          end do
       end if
@@ -277,7 +277,7 @@ SUBROUTINE vol_clu(rho_real,rho_g,s_fac,flag)
 
       e_j = 0.d0
 
-      do ir = 1,nnr
+      do ir = 1,nrxx
    
          v_vol(ir) = 0.d0
 
diff --git a/CPV/wf.f90 b/CPV/wf.f90
index 7440fb43c..e8b2dc914 100644
--- a/CPV/wf.f90
+++ b/CPV/wf.f90
@@ -764,7 +764,7 @@ SUBROUTINE ddyn( m, Omat, Umat, b1, b2, b3 )
   USE wannier_base,     ONLY : wf_friction, nsteps, tolw, adapt, wf_q, &
                                weight, nw, wfdt
   USE cell_base,        ONLY : alat
-  USE constants,        ONLY : tpi, autoaf => BOHR_RADIUS_ANGS
+  USE constants,        ONLY : tpi, bohr_radius_angs
   USE electrons_base,   ONLY : nbsp
   USE control_flags,    ONLY : iprsta
   USE mp_global,        ONLY : me_image
@@ -977,7 +977,7 @@ SUBROUTINE ddyn( m, Omat, Umat, b1, b2, b3 )
   DO i=1, m
      !
      mt=1.D0-DBLE(Oc(:,i,i)*CONJG(Oc(:,i,i)))
-     sp = (alat*autoaf/tpi)**2*SUM(mt*weight)
+     sp = (alat*bohr_radius_angs/tpi)**2*SUM(mt*weight)
      !
      IF(me.EQ.1) THEN
         WRITE(iunit, '(f10.7)') sp
@@ -2492,7 +2492,7 @@ SUBROUTINE wfsteep( m, Omat, Umat, b1, b2, b3 )
   USE wannier_base,           ONLY : nw, weight, nit, tolw, wfdt, maxwfdt, nsd
   USE control_flags,          ONLY : iprsta
   USE cell_base,              ONLY : alat
-  USE constants,              ONLY : tpi, autoaf => BOHR_RADIUS_ANGS
+  USE constants,              ONLY : tpi, bohr_radius_angs
   USE mp_global,              ONLY : me_image
   USE printout_base,          ONLY : printout_base_open, printout_base_unit, &
                                      printout_base_close
@@ -2815,7 +2815,7 @@ SUBROUTINE wfsteep( m, Omat, Umat, b1, b2, b3 )
   DO i=1, m
      !
      mt=1.D0-DBLE(Oc(:,i,i)*CONJG(Oc(:,i,i)))
-     sp = (alat*autoaf/tpi)**2*SUM(mt*weight)
+     sp = (alat*bohr_radius_angs/tpi)**2*SUM(mt*weight)
      !
      IF(me.EQ.1) THEN
         WRITE(iunit, '(f10.7)') sp
@@ -2942,7 +2942,6 @@ SUBROUTINE jacobi_rotation( m, Omat, Umat, b1, b2, b3 )
   !    in the Chair of Atomic Scale Simulation in Lausanne (Switzerland)
   !    under the direction of Prof. Alfredo Pasquarello.
   !
-  REAL(DP), PARAMETER :: autoaf=0.529177d0
   INTEGER, intent(in) :: m
   COMPLEX(DP), DIMENSION(nw, m, m), intent(inout) :: Omat
   REAL(DP), DIMENSION(m, m), intent(inout) :: Umat