[these changes should have been present in a previous commit, but were not]

More removal of old I/O: input_drho, output_drho now use new format. Not sure
they work, though (no test available) and not sure they are used (the previous
version was broken for sure). Fixed a bug added in previous commit.

PW/src/make.depend

@@ -1007,7 +1007,6 @@ io_rho_xml.o : ../../Modules/mp_pools.o
 io_rho_xml.o : ../../Modules/noncol.o
 io_rho_xml.o : ../../Modules/paw_variables.o
 io_rho_xml.o : ../../Modules/recvec.o
-io_rho_xml.o : ../../Modules/xml_io_base.o
 io_rho_xml.o : ../../UtilXlib/mp.o
 io_rho_xml.o : ldaU.o
 io_rho_xml.o : pwcom.o
@@ -1628,11 +1627,17 @@ realus.o : ../../Modules/wavefunctions.o
 realus.o : ../../UtilXlib/mp.o
 realus.o : pwcom.o
 realus.o : scf_mod.o
+remove_atomic_rho.o : ../../Modules/cell_base.o
+remove_atomic_rho.o : ../../Modules/control_flags.o
 remove_atomic_rho.o : ../../Modules/fft_base.o
+remove_atomic_rho.o : ../../Modules/fft_rho.o
+remove_atomic_rho.o : ../../Modules/io_base.o
 remove_atomic_rho.o : ../../Modules/io_files.o
 remove_atomic_rho.o : ../../Modules/io_global.o
 remove_atomic_rho.o : ../../Modules/kind.o
-remove_atomic_rho.o : ../../Modules/xml_io_base.o
+remove_atomic_rho.o : ../../Modules/mp_bands.o
+remove_atomic_rho.o : ../../Modules/mp_pools.o
+remove_atomic_rho.o : ../../Modules/recvec.o
 remove_atomic_rho.o : pwcom.o
 remove_atomic_rho.o : scf_mod.o
 report_mag.o : ../../Modules/constants.o

PW/src/potinit.f90

@@ -83,7 +83,7 @@ SUBROUTINE potinit()
      IF ( .NOT.lforcet ) THEN
         CALL read_scf ( rho, nspin, gamma_only )
         CALL rho_g2r ( dfftp, rho%of_g, rho%of_r )
-        CALL rho_g2r ( dfftp, rho%kin_g, rho%kin_r )
+        IF ( dft_is_meta() ) CALL rho_g2r ( dfftp, rho%kin_g, rho%kin_r )
      ELSE
         !
         ! ... 'force theorem' calculation of MAE: read rho only from previous
@@ -141,13 +141,14 @@ SUBROUTINE potinit()
         IF ( nspin > 1 ) CALL errore &
              ( 'potinit', 'spin polarization not allowed in drho', 1 )
         !
+        filename = TRIM(tmp_dir) // TRIM (prefix) // postfix // input_drho
         CALL read_rhog ( filename, root_bgrp, intra_bgrp_comm, &
              ig_l2g, nspin, v%of_g, gamma_only )
         CALL rho_g2r ( dfftp, v%of_g, v%of_r )
         !
         WRITE( UNIT = stdout, &
                FMT = '(/5X,"a scf correction to at. rho is read from",A)' ) &
-            TRIM( input_drho )
+            TRIM( filename )
         !
         rho%of_r = rho%of_r + v%of_r
         !

PW/src/remove_atomic_rho.f90

@@ -11,33 +11,48 @@ subroutine remove_atomic_rho
   USE io_global, ONLY: stdout
   USE io_files, ONLY: output_drho, tmp_dir, prefix, postfix
   USE kinds, ONLY: DP
+  USE control_flags,    ONLY : gamma_only
+  USE gvect,            ONLY : ig_l2g, mill
   USE fft_base, ONLY: dfftp
+  USE fft_rho, ONLY: rho_r2g
   USE lsda_mod, ONLY: nspin
   USE scf, ONLY: rho
-  USE xml_io_base, ONLY : write_rho
-
+  USE io_base, ONLY : write_rhog
+  USE mp_pools,         ONLY : my_pool_id
+  USE mp_bands,         ONLY : my_bgrp_id, root_bgrp_id, &
+       root_bgrp, intra_bgrp_comm
+  USE cell_base,        ONLY : bg, tpiba
+  !
   implicit none
-  CHARACTER(LEN=256) :: dirname 
+  CHARACTER(LEN=256) :: filename
   real(DP), allocatable :: work (:,:)
   ! workspace, is the difference between the charge density
   ! and the superposition of atomic charges
-
-  allocate ( work( dfftp%nnr, 1 ) )
-  work = 0.d0
+  complex(DP), allocatable :: workc(:,:)
   !
   IF ( nspin > 1 ) CALL errore &
        ( 'remove_atomic_rho', 'spin polarization not allowed in drho', 1 )
 
   WRITE( stdout, '(/5x,"remove atomic charge density from scf rho")')
   !
-  !     subtract the old atomic charge density
+  !     subtract the old atomic charge density (FIXME: in real space)
   !
+  allocate ( work( dfftp%nnr, 1 ) )
+  work = 0.d0
   call atomic_rho (work, nspin)
-  !
   work = rho%of_r - work
   !
-  dirname = TRIM(tmp_dir) // TRIM(prefix) // postfix
-  call write_rho ( dirname, work, 1, output_drho )
+  !      FIXME: move to G-space
+  !
+  allocate ( workc( dfftp%ngm, 1) )
+  CALL rho_r2g ( dfftp, work, workc )
+  deallocate(work)
+  !
+  filename = TRIM(tmp_dir) // TRIM(prefix) // postfix // output_drho
+  IF ( my_pool_id == 0 .AND. my_bgrp_id == root_bgrp_id ) &
+       CALL write_rhog( filename, root_bgrp, intra_bgrp_comm, &
+       bg(:,1)*tpiba, bg(:,2)*tpiba, bg(:,3)*tpiba, &
+       gamma_only, mill, ig_l2g, workc )
   !
   deallocate(work)
   return