Modification by M.Cococcioni

Here is the accompanying explanation:
"What I achieved is that at the last step of a bfgs calculation, when it makes a final scf restarting from scratch (and latest atomic positions and cell) it reads the occupations from file. This way it does not loose the character of the ground state it had achieved during the DFT+U(+V) scf and vc-relax. Sometimes it is qualitatively different from a DFT one (e.g. it has a different order of states) and washing this information is not good. In other words, the Hamiltonian and its ground state are functions of also the Hubbard parameters so one need to preserve the ground state they achieved. removing this information is like resetting the atomic position to the initial values after each vc-relax..."

Doc/release-notes

@@ -1,6 +1,9 @@
 New in dev version:
   * Support for CMake (F. Ficarelli and D. Cesarini, CINECA, with help from
     Ye Luo, P. Delugas, S. Gsaenger)
+  * In vc-relax with Hubbard corrections, the final SCF calculation is done by 
+    reading atomic occupations from file produced during the vc-relax 
+    (rather then recomputing them from scratch).
 
 Fixed in dev version:
   * Some linkers yield "missing references to ddot_" in libbeef

PW/src/run_pwscf.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2013-2017 Quantum ESPRESSO group
+! Copyright (C) 2013-2020 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,
@@ -33,6 +33,8 @@ SUBROUTINE run_pwscf( exit_status )
   !! @endnote
   !!
   !
+  USE kinds,                ONLY : DP
+  USE mp,                   ONLY : mp_bcast, mp_sum
   USE io_global,            ONLY : stdout, ionode, ionode_id
   USE parameters,           ONLY : ntypx, npk
   USE upf_params,           ONLY : lmaxx
@@ -53,6 +55,7 @@ SUBROUTINE run_pwscf( exit_status )
   USE qexsd_module,         ONLY : qexsd_set_status
   USE funct,                ONLY : dft_is_hybrid, stop_exx 
   USE beef,                 ONLY : beef_energies
+  USE ldaU,                 ONLY : lda_plus_u
   !
   IMPLICIT NONE
   !
@@ -228,8 +231,13 @@ SUBROUTINE run_pwscf( exit_status )
            !
            lbfgs=.FALSE.; lmd=.FALSE.
            WRITE( UNIT = stdout, FMT=9020 ) 
+           !
            CALL reset_gvectors( )
            !
+           ! ... read atomic occupations for DFT+U(+V)
+           !
+           IF ( lda_plus_u ) CALL read_ns()
+           !
         ELSE IF ( ions_status == 2 ) THEN
            !
            ! ... check whether nonzero magnetization is real

PW/src/write_ns.f90

@@ -671,3 +671,94 @@ SUBROUTINE write_nsg
   !
 END SUBROUTINE write_nsg 
 
+!-----------------------------------------------------------------------
+SUBROUTINE read_ns()
+  !---------------------------------------------------------------------
+  !
+  ! This routine was written for the final SCF after vc-relax (M. Cococcioni).
+  ! The occupations ns/nsg also need to be read in order to reproduce
+  ! the right electronic ground state. When using Hubbard corrections
+  ! this might be different (i.e. have different ordering of states) 
+  ! from that simply obtained from the superposition of free ions.
+  ! In other words the KS Hamiltonian (and its ground state) is also 
+  ! functional of the Hubbard interaction parameters.
+  !
+  USE kinds,              ONLY : DP
+  USE mp,                 ONLY : mp_bcast
+  USE mp_images,          ONLY : intra_image_comm
+  USE io_global,          ONLY : ionode, ionode_id
+  USE scf,                ONLY : rho, v
+  USE ldaU,               ONLY : lda_plus_u_kind, nsg, v_nsg, hub_back
+  USE noncollin_module,   ONLY : noncolin
+  USE io_files,           ONLY : restart_dir
+  !
+  IMPLICIT NONE
+  INTEGER :: iunocc, iunocc1, ierr
+  CHARACTER (LEN=256) :: dirname
+  REAL(DP) :: eth, eth1
+  !
+  dirname = restart_dir()
+  !
+  IF ( ionode ) THEN
+     !
+     OPEN ( NEWUNIT=iunocc, FILE = TRIM(dirname) // 'occup.txt', &
+            FORM='formatted', STATUS='old', IOSTAT=ierr )
+     IF (lda_plus_u_kind.EQ.0) THEN
+        READ( UNIT = iunocc, FMT = *, iostat = ierr ) rho%ns
+        IF (hub_back) THEN
+           READ( UNIT = iunocc, FMT = * , iostat = ierr) rho%nsb
+        ENDIF
+     ELSEIF (lda_plus_u_kind.EQ.1) THEN
+        IF (noncolin) THEN
+           READ( UNIT = iunocc, FMT = *, iostat = ierr ) rho%ns_nc
+        ELSE
+           READ( UNIT = iunocc, FMT = *, iostat = ierr ) rho%ns
+        ENDIF
+     ELSEIF (lda_plus_u_kind.EQ.2) THEN
+        READ( UNIT = iunocc, FMT = * , iostat = ierr) nsg
+     ENDIF
+     CLOSE(UNIT=iunocc,STATUS='keep')
+     !
+  ELSE
+     !
+     IF (lda_plus_u_kind.EQ.0) THEN
+        rho%ns(:,:,:,:) = 0.D0
+        IF (hub_back) rho%nsb(:,:,:,:) = 0.D0
+     ELSEIF (lda_plus_u_kind.EQ.1) THEN
+        IF (noncolin) THEN
+           rho%ns_nc(:,:,:,:) = 0.D0
+        ELSE
+           rho%ns(:,:,:,:) = 0.D0
+        ENDIF
+     ELSEIF (lda_plus_u_kind.EQ.2) THEN
+        nsg(:,:,:,:,:) = (0.d0, 0.d0)
+     ENDIF
+     !
+  ENDIF
+  !
+  CALL mp_bcast( ierr, ionode_id, intra_image_comm )
+  !
+  IF (lda_plus_u_kind.EQ.0) THEN
+     CALL mp_bcast(rho%ns, ionode_id, intra_image_comm)
+     CALL v_hubbard (rho%ns, v%ns, eth)
+     IF (hub_back) THEN
+        CALL mp_bcast(rho%nsb, ionode_id, intra_image_comm)
+        CALL v_hubbard_b (rho%nsb, v%nsb, eth1)
+        eth = eth + eth1
+     ENDIF
+  ELSEIF (lda_plus_u_kind.EQ.1) THEN
+     IF (noncolin) THEN
+        CALL mp_bcast(rho%ns_nc, ionode_id, intra_image_comm)
+        CALL v_hubbard_full_nc (rho%ns_nc, v%ns_nc, eth)
+     ELSE
+        CALL mp_bcast(rho%ns, ionode_id, intra_image_comm)
+        CALL v_hubbard_full (rho%ns, v%ns, eth)
+     ENDIF
+  ELSEIF (lda_plus_u_kind.EQ.2) THEN
+     CALL mp_bcast(nsg, ionode_id, intra_image_comm)
+     CALL v_hubbard_extended (nsg, v_nsg, eth)
+  ENDIF
+  !
+  RETURN
+  !
+END SUBROUTINE read_ns