The charge density is written and read to the XML restart file

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@2851 c92efa57-630b-4861-b058-cf58834340f0

PW/pw_restart.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2005 Quantum-ESPRESSO group
+! Copyright (C) 2005-2006 Quantum-ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
@@ -97,7 +97,7 @@ MODULE pw_restart
       CHARACTER(LEN=256)    :: dirname, filename, file_pseudo, rho_file_base
       CHARACTER(LEN=80)     :: bravais_lattice
       INTEGER               :: i, ig, ik, ngg, ierr, ipol, ik_eff, num_k_points
-      INTEGER,  ALLOCATABLE :: kisort(:)
+      !!!INTEGER,  ALLOCATABLE :: kisort(:)
       INTEGER               :: npool, nkbl, nkl, nkr, npwx_g
       INTEGER               :: ike, iks, npw_g, ispin
       INTEGER,  ALLOCATABLE :: ngk_g(:)
@@ -204,23 +204,23 @@ MODULE pw_restart
       !
       ! ... build the G+k array indexes
       !
-      ALLOCATE( kisort( npwx ) )
+      !!!ALLOCATE( kisort( npwx ) )
       !
-      DO ik = 1, nks
+      !!!DO ik = 1, nks
          !
-         kisort = 0
-         npw    = npwx
+         !!!kisort = 0
+         !!!npw    = npwx
          !
-         CALL gk_sort( xk(1,ik+iks-1), ngm, g, &
-                       ecutwfc/tpiba2, npw, kisort(1), g2kin )
+         !!!CALL gk_sort( xk(1,ik+iks-1), ngm, g, &
+         !!!              ecutwfc/tpiba2, npw, kisort(1), g2kin )
          !
-         CALL gk_l2gmap( ngm, ig_l2g(1), npw, kisort(1), igk_l2g(1,ik) )
+         !!!CALL gk_l2gmap( ngm, ig_l2g(1), npw, kisort(1), igk_l2g(1,ik) )
          !
-         ngk(ik) = npw
+         !!!ngk(ik) = npw
          !
-      END DO
+      !!!END DO
       !
-      DEALLOCATE( kisort )
+      !!!DEALLOCATE( kisort )
       !
       ! ... compute the global number of G+k vectors for each k point
       !
@@ -398,6 +398,36 @@ MODULE pw_restart
          CALL write_rho_xml( rho_file_base, rho(:,2), &
                              nr1, nr2, nr3, nrx1, nrx2, dfftp%ipp, dfftp%npp )
          !
+      ELSE IF ( nspin == 4 ) THEN
+         !
+         CALL write_rho_xml( rho_file_base, rho(:,1), nr1, &
+                             nr2, nr3, nrx1, nrx2, dfftp%ipp, dfftp%npp )
+         !
+         rho_file_base = 'magnetization.x'
+         IF ( ionode ) &
+            CALL iotk_link( iunpun, "MAG_X", rho_file_base, &
+                            CREATE = .FALSE., BINARY = .FALSE. )
+         rho_file_base = TRIM( dirname ) // '/' // TRIM( rho_file_base )
+         CALL write_rho_xml( rho_file_base, rho(:,2), &
+                             nr1, nr2, nr3, nrx1, nrx2, dfftp%ipp, dfftp%npp )
+         !
+         rho_file_base = 'magnetization.y'
+         IF ( ionode ) &
+            CALL iotk_link( iunpun, "MAG_Y", rho_file_base, &
+                            CREATE = .FALSE., BINARY = .FALSE. )
+         rho_file_base = TRIM( dirname ) // '/' // TRIM( rho_file_base )
+         CALL write_rho_xml( rho_file_base, rho(:,3), &
+                             nr1, nr2, nr3, nrx1, nrx2, dfftp%ipp, dfftp%npp )
+         !
+         !
+         rho_file_base = 'magnetization.z'
+         IF ( ionode ) &
+            CALL iotk_link( iunpun, "MAG_Z", rho_file_base, &
+                            CREATE = .FALSE., BINARY = .FALSE. )
+         rho_file_base = TRIM( dirname ) // '/' // TRIM( rho_file_base )
+         CALL write_rho_xml( rho_file_base, rho(:,4), &
+                             nr1, nr2, nr3, nrx1, nrx2, dfftp%ipp, dfftp%npp )
+         !
       END IF
       !
       IF ( ionode ) THEN
@@ -725,7 +755,7 @@ MODULE pw_restart
       CHARACTER(LEN=256) :: dirname      
       LOGICAL            :: lexist, lcell, lpw, lions, lspin, &
                             lxc, locc, lbz, lbs, lwfc, lgvec, &
-                            lsymm, lph
+                            lsymm, lph, lrho
       !
       !
       ierr = 0
@@ -761,6 +791,7 @@ MODULE pw_restart
       lgvec = .FALSE.
       lsymm = .FALSE.
       lph   = .FALSE.
+      lrho  = .FALSE.
       !
       SELECT CASE( what )
       CASE( 'dim' )
@@ -778,6 +809,10 @@ MODULE pw_restart
          lcell = .TRUE.
          lions = .TRUE.
          !
+      CASE( 'rho' )
+         !
+         lrho  = .TRUE.
+         !
       CASE( 'wave' )
          !
          lpw   = .TRUE.
@@ -811,6 +846,7 @@ MODULE pw_restart
          lgvec = .TRUE.
          lsymm = .TRUE.
          lph   = .TRUE.
+         lrho  = .TRUE.
          !
       CASE( 'reset' )
          !
@@ -895,6 +931,13 @@ MODULE pw_restart
          !
       END IF
       !
+      IF ( lrho ) THEN
+         !
+         CALL read_rho( dirname, ierr )
+         IF ( ierr > 0 ) RETURN
+         !
+      END IF
+      !
       RETURN
       !
     END SUBROUTINE pw_readfile
@@ -2192,6 +2235,61 @@ MODULE pw_restart
     END SUBROUTINE read_phonon
     !
     !------------------------------------------------------------------------
+    SUBROUTINE read_rho( dirname, ierr )
+      !------------------------------------------------------------------------
+      !
+      USE gvect,                ONLY : nr1, nr2, nr3, nrx1, nrx2
+      USE lsda_mod,             ONLY : nspin
+      USE scf,                  ONLY : rho
+      USE sticks,               ONLY : dfftp
+      !
+      IMPLICIT NONE
+      !
+      CHARACTER(LEN=*), INTENT(IN)  :: dirname
+      INTEGER,          INTENT(OUT) :: ierr
+      !
+      CHARACTER(LEN=256) :: rho_file_base
+      !
+      IF (nspin == 1) THEN
+         !
+         rho_file_base = TRIM( dirname ) // '/charge-density'
+         CALL read_rho_xml( rho_file_base, rho(:,1), &
+              nr1, nr2, nr3, nrx1, nrx2,  dfftp%ipp, dfftp%npp )
+         !
+      ELSE IF (nspin == 2) THEN
+         !
+         rho_file_base = TRIM( dirname ) // '/charge-density-up'
+         CALL read_rho_xml( rho_file_base, rho(:,1), &
+              nr1, nr2, nr3, nrx1, nrx2,  dfftp%ipp, dfftp%npp )
+         !
+         rho_file_base = TRIM( dirname ) // '/charge-density-dw'
+         CALL read_rho_xml( rho_file_base, rho(:,2), &
+              nr1, nr2, nr3, nrx1, nrx2,  dfftp%ipp, dfftp%npp )
+         !
+      ELSE IF (nspin == 4) THEN
+         !
+         rho_file_base = TRIM( dirname ) // '/charge-density'
+         CALL read_rho_xml( rho_file_base, rho(:,1), &
+              nr1, nr2, nr3, nrx1, nrx2,  dfftp%ipp, dfftp%npp )
+         !
+         rho_file_base = TRIM( dirname ) // '/magnetization.x'
+         CALL read_rho_xml( rho_file_base, rho(:,2), &
+              nr1, nr2, nr3, nrx1, nrx2,  dfftp%ipp, dfftp%npp )
+         !
+         rho_file_base = TRIM( dirname ) // '/magnetization.y'
+         CALL read_rho_xml( rho_file_base, rho(:,3), &
+              nr1, nr2, nr3, nrx1, nrx2,  dfftp%ipp, dfftp%npp )
+         !
+         rho_file_base = TRIM( dirname ) // '/magnetization.z'
+         CALL read_rho_xml( rho_file_base, rho(:,4), &
+              nr1, nr2, nr3, nrx1, nrx2,  dfftp%ipp, dfftp%npp )
+         !
+      END IF
+      ierr = 0
+      !
+    END SUBROUTINE read_rho
+    !
+    !------------------------------------------------------------------------
     SUBROUTINE read_( dirname, ierr )
       !------------------------------------------------------------------------
       !

PW/read_file.f90

@@ -17,7 +17,7 @@ SUBROUTINE read_file()
   USE kinds,            ONLY : DP
   USE ions_base,        ONLY : nat, nsp, ityp, tau, if_pos
   USE basis,            ONLY : natomwfc
-  USE cell_base,        ONLY : tpiba2, at, bg
+  USE cell_base,        ONLY : tpiba2, alat,omega,  at, bg
   USE force_mod,        ONLY : force
   USE klist,            ONLY : nkstot, nks, xk, wk
   USE lsda_mod,         ONLY : lsda, nspin, current_spin, isk
@@ -26,11 +26,12 @@ SUBROUTINE read_file()
   USE ktetra,           ONLY : tetra, ntetra 
   USE extfield,         ONLY : forcefield, tefield
   USE cellmd,           ONLY : cell_factor, lmovecell
-  USE gvect,            ONLY : gg, ecutwfc, ngm, g, nr1, nr2, nr3, &
-                               eigts1, eigts2, eigts3
+  USE gvect,            ONLY : gg, ecutwfc, ngm, g, nr1, nr2, nr3, nrxx,&
+                               nrx1, nrx2, nrx3, eigts1, eigts2, eigts3, &
+                               nl, gstart
   USE gsmooth,          ONLY : ngms, nls, nrx1s, nr1s, nr2s, nr3s
   USE spin_orb,         ONLY : so
-  USE scf,              ONLY : rho, vr
+  USE scf,              ONLY : rho, rho_core, vr
   USE vlocal,           ONLY : strf
   USE io_files,         ONLY : tmp_dir, prefix, iunpun
   USE restart_module,   ONLY : readfile_new
@@ -43,7 +44,7 @@ SUBROUTINE read_file()
   IMPLICIT NONE
   !
   INTEGER               :: i, ik, ibnd, nb, nt, ios, ierr
-  REAL(DP)              :: rdum(1,1)
+  REAL(DP)              :: rdum(1,1), ehart, etxc, vtxc, etotefield, charge
   !
   !
   !
@@ -180,6 +181,29 @@ SUBROUTINE read_file()
   CALL allocate_nlpot()
   CALL allocate_wfc()
   !
+  ! ... read the charge density
+  !
+  CALL pw_readfile( 'rho', ierr )
+  !
+  ! ... re-calculate the local part of the pseudopotential vltot
+  ! ... and the core correction charge (if any) - This is done here
+  ! ... for compatibility with the previous version of read_file
+  !
+  CALL init_vloc()
+  !
+  CALL struc_fact( nat, tau, nsp, ityp, ngm, g, bg, &
+                   nr1, nr2, nr3, strf, eigts1, eigts2, eigts3 )
+  !
+  CALL setlocal()
+  !
+  CALL set_rhoc()
+  !
+  ! ... recalculate the potential
+  !
+  CALL v_of_rho( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3,   &
+       nrxx, nl, ngm, gstart, nspin, g, gg, alat, omega, &
+       ehart, etxc, vtxc, etotefield, charge, vr )
+  !
 #else
   !
   !-------------------------------------------------------------------------------
@@ -258,8 +282,6 @@ SUBROUTINE read_file()
   CALL poolscatter( nbnd , nkstot, et, nks, et )
   CALL poolscatter( nbnd , nkstot, wg, nks, wg )
   !
-#endif
-  !
   ! ... read the charge density
   !
   CALL io_pot( - 1, 'rho', rho, nspin )
@@ -281,6 +303,8 @@ SUBROUTINE read_file()
   !
   CALL set_rhoc()
   !
+#endif
+  !
   RETURN
   !
   CONTAINS