diff --git a/CPV/src/chargedensity.f90 b/CPV/src/chargedensity.f90
index a19941a79..e8d5a911b 100644
--- a/CPV/src/chargedensity.f90
+++ b/CPV/src/chargedensity.f90
@@ -195,7 +195,7 @@
                  TRIM(tmp_dir), TRIM(prefix), ndr
             CALL read_rhog ( dirname, root_bgrp, intra_bgrp_comm, &
                  ig_l2g, nspin, rhog )
-            CALL rho_g2r ( rhog, rhor )
+            CALL rho_g2r ( dfftp, rhog, rhor )
 #endif
             rhopr = rhor
             first = .FALSE.
@@ -255,7 +255,7 @@
          !
          rhog(dffts%ngm+1:,:) = 0.0d0
          !
-         CALL rho_g2r( rhog, rhor )
+         CALL rho_g2r( dfftp, rhog, rhor )
          !
          IF ( dft_is_meta() ) THEN
             CALL kedtauofr_meta( c_bgrp ) ! METAGGA
diff --git a/CPV/src/cplib_meta.f90 b/CPV/src/cplib_meta.f90
index 869481f98..4a381a351 100644
--- a/CPV/src/cplib_meta.f90
+++ b/CPV/src/cplib_meta.f90
@@ -192,7 +192,7 @@
 !
       kedtaug(dffts%ngm+1:,:) = 0.0d0
 
-      CALL rho_g2r( kedtaug, kedtaur )
+      CALL rho_g2r( dfftp, kedtaug, kedtaur )
 !
       return
     end subroutine kedtauofr_meta
@@ -258,7 +258,7 @@
 !
       CALL rho_r2g( dfftp, kedtaur, kedtaug )
 !
-      CALL smooth_rho_g2r ( kedtaug, kedtaus )
+      CALL rho_g2r ( dffts, kedtaug, kedtaus )
 
       !calculate dkedxc in real space on smooth grids  !metagga
       if(tpre) then
diff --git a/CPV/src/vofrho.f90 b/CPV/src/vofrho.f90
index b085c966a..9fa81cb8f 100644
--- a/CPV/src/vofrho.f90
+++ b/CPV/src/vofrho.f90
@@ -474,7 +474,7 @@ SUBROUTINE vofrho_x( nfi, rhor, drhor, rhog, drhog, rhos, rhoc, tfirst, &
 !     fourier transform of total potential to r-space (dense grid)
 !     -------------------------------------------------------------------
       
-      CALL rho_g2r( rhog, rhor )
+      CALL rho_g2r( dfftp, rhog, rhor )
 
       IF(nspin.EQ.1) THEN
          vave=SUM(rhor(:,1))/DBLE( dfftp%nr1* dfftp%nr2* dfftp%nr3)
@@ -489,7 +489,7 @@ SUBROUTINE vofrho_x( nfi, rhor, drhor, rhog, drhog, rhos, rhoc, tfirst, &
       !
       !     fourier transform of total potential to r-space (smooth grid)
       !
-      CALL smooth_rho_g2r ( rhog, rhos )
+      CALL rho_g2r ( dffts, rhog, rhos )
 
       IF( dft_is_meta() ) CALL vofrho_meta( )
 
diff --git a/CPV/src/vol_clu.f90 b/CPV/src/vol_clu.f90
index 0c602401e..19801ba80 100644
--- a/CPV/src/vol_clu.f90
+++ b/CPV/src/vol_clu.f90
@@ -244,7 +244,7 @@ SUBROUTINE vol_clu(rho_real,rho_g,flag)
       if (abisur)                                                       &
      &   call gradrho(nspin,rhotmp,drho,d2rho,dxdyrho,dxdzrho,dydzrho)
 
-      CALL rho_g2r( rhotmp, rho_gaus )
+      CALL rho_g2r( dfftp, rhotmp, rho_gaus )
       deallocate(rhotmp)
 
       e_j = 0.d0
diff --git a/Modules/fft_rho.f90 b/Modules/fft_rho.f90
index 7b0faec41..d771499ef 100644
--- a/Modules/fft_rho.f90
+++ b/Modules/fft_rho.f90
@@ -17,7 +17,7 @@ MODULE fft_rho
   !
   IMPLICIT NONE
   PRIVATE
-  PUBLIC :: rho_r2g, rho_g2r, smooth_rho_g2r
+  PUBLIC :: rho_r2g, rho_g2r
   !
   INTERFACE rho_g2r
     MODULE PROCEDURE rho_g2r_x, rho_g2r_sum_components
@@ -75,10 +75,11 @@ CONTAINS
 
   END SUBROUTINE rho_r2g
   !
-  SUBROUTINE rho_g2r_x ( rhog, rhor )
-    USE fft_base,       ONLY: dfftp
+  SUBROUTINE rho_g2r_x ( desc, rhog, rhor )
+    USE fft_types,              ONLY: fft_type_descriptor
     USE fft_helper_subroutines, ONLY: fftx_threed2oned, fftx_oned2threed
     !
+    TYPE(fft_type_descriptor), INTENT(in) :: desc
     COMPLEX(dp), INTENT(in ):: rhog(:,:)
     REAL(dp),    INTENT(out):: rhor(:,:)
     !
@@ -89,24 +90,24 @@ CONTAINS
 
     nspin= SIZE (rhog, 2)
 
-    ALLOCATE( psi( dfftp%nnr ) )
+    ALLOCATE( psi( desc%nnr ) )
     IF ( gamma_only ) THEN
        IF( nspin == 1 ) THEN
           iss=1
-          CALL fftx_oned2threed( dfftp, psi, rhog(:,iss) )
-          CALL invfft('Rho',psi, dfftp )
+          CALL fftx_oned2threed( desc, psi, rhog(:,iss) )
+          CALL invfft('Rho',psi, desc )
 !$omp parallel do
-          DO ir=1,dfftp%nnr
+          DO ir=1,desc%nnr
              rhor(ir,iss)=DBLE(psi(ir))
           END DO
 !$omp end parallel do
        ELSE
           isup=1
           isdw=2
-          CALL fftx_oned2threed( dfftp, psi, rhog(:,isup), rhog(:,isdw) )
-          CALL invfft('Rho',psi, dfftp )
+          CALL fftx_oned2threed( desc, psi, rhog(:,isup), rhog(:,isdw) )
+          CALL invfft('Rho',psi, desc )
 !$omp parallel do
-          DO ir=1,dfftp%nnr
+          DO ir=1,desc%nnr
              rhor(ir,isup)= DBLE(psi(ir))
              rhor(ir,isdw)=AIMAG(psi(ir))
           END DO
@@ -116,10 +117,10 @@ CONTAINS
     ELSE
        !
        DO iss=1, nspin
-          CALL fftx_oned2threed( dfftp, psi, rhog(:,iss) )
-          CALL invfft('Rho',psi, dfftp )
+          CALL fftx_oned2threed( desc, psi, rhog(:,iss) )
+          CALL invfft('Rho',psi, desc )
 !$omp parallel do
-          DO ir=1,dfftp%nnr
+          DO ir=1,desc%nnr
              rhor(ir,iss)=DBLE(psi(ir))
           END DO
 !$omp end parallel do
@@ -130,10 +131,11 @@ CONTAINS
 
   END SUBROUTINE rho_g2r_x
   !
-  SUBROUTINE rho_g2r_sum_components ( rhog, rhor )
-    USE fft_base,       ONLY: dfftp
+  SUBROUTINE rho_g2r_sum_components ( desc, rhog, rhor )
+    USE fft_types,              ONLY: fft_type_descriptor
     USE fft_helper_subroutines, ONLY: fftx_threed2oned, fftx_oned2threed
     !
+    TYPE(fft_type_descriptor), INTENT(in) :: desc
     COMPLEX(dp), INTENT(in ):: rhog(:,:)
     REAL(dp),    INTENT(out):: rhor(:)
     !
@@ -144,24 +146,24 @@ CONTAINS
 
     nspin= SIZE (rhog, 2)
 
-    ALLOCATE( psi( dfftp%nnr ) )
+    ALLOCATE( psi( desc%nnr ) )
     IF ( gamma_only ) THEN
        IF( nspin == 1 ) THEN
           iss=1
-          CALL fftx_oned2threed( dfftp, psi, rhog(:,iss) )
-          CALL invfft('Rho',psi, dfftp )
+          CALL fftx_oned2threed( desc, psi, rhog(:,iss) )
+          CALL invfft('Rho',psi, desc )
 !$omp parallel do
-          DO ir=1,dfftp%nnr
+          DO ir=1,desc%nnr
              rhor(ir)=DBLE(psi(ir))
           END DO
 !$omp end parallel do
        ELSE
           isup=1
           isdw=2
-          CALL fftx_oned2threed( dfftp, psi, rhog(:,isup), rhog(:,isdw) )
-          CALL invfft('Rho',psi, dfftp )
+          CALL fftx_oned2threed( desc, psi, rhog(:,isup), rhog(:,isdw) )
+          CALL invfft('Rho',psi, desc )
 !$omp parallel do
-          DO ir=1,dfftp%nnr
+          DO ir=1,desc%nnr
              rhor(ir)= DBLE(psi(ir))+AIMAG(psi(ir))
           END DO
 !$omp end parallel do
@@ -170,17 +172,17 @@ CONTAINS
     ELSE
        !
        DO iss=1, nspin
-          CALL fftx_oned2threed( dfftp, psi, rhog(:,iss) )
-          CALL invfft('Rho',psi, dfftp )
+          CALL fftx_oned2threed( desc, psi, rhog(:,iss) )
+          CALL invfft('Rho',psi, desc )
           IF( iss == 1 ) THEN
 !$omp parallel do
-             DO ir=1,dfftp%nnr
+             DO ir=1,desc%nnr
                 rhor(ir)=DBLE(psi(ir))
              END DO
 !$omp end parallel do
           ELSE
 !$omp parallel do
-             DO ir=1,dfftp%nnr
+             DO ir=1,desc%nnr
                 rhor(ir)=rhor(ir) + DBLE(psi(ir))
              END DO
 !$omp end parallel do
@@ -192,60 +194,4 @@ CONTAINS
 
   END SUBROUTINE rho_g2r_sum_components
 
-  SUBROUTINE smooth_rho_g2r ( rhog, rhor )
-    USE fft_base,       ONLY: dffts
-    USE fft_helper_subroutines, ONLY: fftx_threed2oned, fftx_oned2threed
-    !
-    COMPLEX(dp), INTENT(in ):: rhog(:,:)
-    REAL(dp),    INTENT(out):: rhor(:,:)
-    !
-    INTEGER :: ir, ig, iss, isup, isdw
-    INTEGER :: nspin
-    COMPLEX(dp), PARAMETER :: ci=(0.0_dp, 1.0_dp)
-    COMPLEX(dp), ALLOCATABLE :: psi(:)
-
-    nspin= SIZE (rhog, 2)
-
-    ALLOCATE( psi( dffts%nnr ) )
-    IF ( gamma_only ) THEN
-       IF( nspin == 1 ) THEN
-          iss=1
-          CALL fftx_oned2threed( dffts, psi, rhog(:,iss) )
-          CALL invfft('Rho',psi, dffts )
-!$omp parallel do
-          DO ir=1,dffts%nnr
-             rhor(ir,iss)=DBLE(psi(ir))
-          END DO
-!$omp end parallel do
-       ELSE
-          isup=1
-          isdw=2
-          CALL fftx_oned2threed( dffts, psi, rhog(:,isup), rhog(:,isdw) )
-          CALL invfft('Rho',psi, dffts )
-!$omp parallel do
-          DO ir=1,dffts%nnr
-             rhor(ir,isup)= DBLE(psi(ir))
-             rhor(ir,isdw)=AIMAG(psi(ir))
-          END DO
-!$omp end parallel do
-       ENDIF
-       !
-    ELSE
-       !
-       DO iss=1, nspin
-          CALL fftx_oned2threed( dffts, psi, rhog(:,iss) )
-          CALL invfft('Rho',psi, dffts )
-!$omp parallel do
-          DO ir=1,dffts%nnr
-             rhor(ir,iss)=DBLE(psi(ir))
-          END DO
-!$omp end parallel do
-       END DO
-    END IF
-    
-    DEALLOCATE( psi )
-
-  END SUBROUTINE smooth_rho_g2r
-
-
 END MODULE fft_rho
diff --git a/PW/src/potinit.f90 b/PW/src/potinit.f90
index eb42732ed..fdeabd5b0 100644
--- a/PW/src/potinit.f90
+++ b/PW/src/potinit.f90
@@ -89,7 +89,7 @@ SUBROUTINE potinit()
      IF ( .NOT.lforcet ) THEN
         CALL read_scf ( rho, nspin, gamma_only )
 #if !defined (__OLDXML)
-        CALL rho_g2r ( rho%of_g, rho%of_r )
+        CALL rho_g2r ( dfftp, rho%of_g, rho%of_r )
 #endif
      ELSE
         !
@@ -101,7 +101,7 @@ SUBROUTINE potinit()
 #else
         CALL read_rhog ( dirname, root_bgrp, intra_bgrp_comm, &
              ig_l2g, nspin, rho%of_g, gamma_only )
-        CALL rho_g2r ( rho%of_g, rho%of_r )
+        CALL rho_g2r ( dfftp, rho%of_g, rho%of_r )
 #endif
         CALL nc_magnetization_from_lsda ( dfftp%nnr, nspin, rho%of_r )
      END IF
@@ -157,7 +157,7 @@ SUBROUTINE potinit()
 #else
         CALL read_rhog ( dirname, root_bgrp, intra_bgrp_comm, &
              ig_l2g, nspin, v%of_g, gamma_only )
-        CALL rho_g2r ( v%of_g, v%of_r )
+        CALL rho_g2r ( dfftp, v%of_g, v%of_r )
 #endif
         !
         WRITE( UNIT = stdout, &
diff --git a/PW/src/read_file_new.f90 b/PW/src/read_file_new.f90
index 663718355..f60514e0f 100644
--- a/PW/src/read_file_new.f90
+++ b/PW/src/read_file_new.f90
@@ -305,7 +305,7 @@ SUBROUTINE read_xml_file ( )
 #if ! defined (__OLDXML)
   ! FIXME: for compatibility. rho was previously read and written in real space
   ! FIXME: now it is in G space - to be removed together with old format
-  CALL rho_g2r ( rho%of_g, rho%of_r )
+  CALL rho_g2r ( dfftp, rho%of_g, rho%of_r )
 #endif
   !
   ! ... re-calculate the local part of the pseudopotential vltot