mirror of https://gitlab.com/QEF/q-e.git
Some simplifications in the logic of "electron":
- non-scf calculation done in a separate routine - check on insufficient diagonalization threshold less obscure git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1889 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
parent
9ee776b63e
commit
818ee4b111
242
PW/electrons.f90
242
PW/electrons.f90
|
@ -57,7 +57,6 @@ SUBROUTINE electrons()
|
||||||
niter_with_fixed_ns, Hubbard_lmax, &
|
niter_with_fixed_ns, Hubbard_lmax, &
|
||||||
lda_plus_u
|
lda_plus_u
|
||||||
USE extfield, ONLY : tefield, etotefield
|
USE extfield, ONLY : tefield, etotefield
|
||||||
USE bp, ONLY : lberry
|
|
||||||
USE wavefunctions_module, ONLY : evc, evc_nc
|
USE wavefunctions_module, ONLY : evc, evc_nc
|
||||||
USE noncollin_module, ONLY : noncolin, npol, magtot_nc
|
USE noncollin_module, ONLY : noncolin, npol, magtot_nc
|
||||||
USE noncollin_module, ONLY : factlist, pointlist, pointnum, mcons,&
|
USE noncollin_module, ONLY : factlist, pointlist, pointnum, mcons,&
|
||||||
|
@ -101,13 +100,12 @@ SUBROUTINE electrons()
|
||||||
etotefield_new, &!
|
etotefield_new, &!
|
||||||
charge_new !
|
charge_new !
|
||||||
REAL (KIND=DP) :: &
|
REAL (KIND=DP) :: &
|
||||||
ethr_min ! minimal threshold for diagonalization at the first scf
|
tr2_min ! estimated error on the energy coming from diagonalization
|
||||||
! iteration
|
|
||||||
REAL (KIND=DP), ALLOCATABLE :: &
|
REAL (KIND=DP), ALLOCATABLE :: &
|
||||||
wg_g(:,:) ! temporary array used to recover from pools array wg,
|
wg_g(:,:) ! temporary array used to recover from pools array wg,
|
||||||
! and then print occupations on stdout
|
! and then print occupations on stdout
|
||||||
LOGICAL :: &
|
LOGICAL :: &
|
||||||
exst
|
exst, first
|
||||||
REAL (KIND=DP) :: dr2old, lambda0
|
REAL (KIND=DP) :: dr2old, lambda0
|
||||||
!
|
!
|
||||||
! ... external functions
|
! ... external functions
|
||||||
|
@ -152,69 +150,7 @@ SUBROUTINE electrons()
|
||||||
!
|
!
|
||||||
IF ( .NOT. lscf ) THEN
|
IF ( .NOT. lscf ) THEN
|
||||||
!
|
!
|
||||||
WRITE( stdout, 9002 )
|
CALL non_scf()
|
||||||
!
|
|
||||||
#if defined (FLUSH)
|
|
||||||
CALL flush( stdout )
|
|
||||||
#endif
|
|
||||||
!
|
|
||||||
IF ( imix >= 0 ) rho_save = rho
|
|
||||||
!
|
|
||||||
iter = iter + 1
|
|
||||||
!
|
|
||||||
! ... diagonalization of the KS hamiltonian
|
|
||||||
!
|
|
||||||
CALL c_bands( iter, ik_, dr2 )
|
|
||||||
!
|
|
||||||
conv_elec = .TRUE.
|
|
||||||
!
|
|
||||||
CALL poolrecover( et, nbnd, nkstot, nks )
|
|
||||||
!
|
|
||||||
tcpu = get_clock( 'PWSCF' )
|
|
||||||
WRITE( stdout, 9000 ) tcpu
|
|
||||||
!
|
|
||||||
WRITE( stdout, 9102 )
|
|
||||||
!
|
|
||||||
! ... write band eigenvalues
|
|
||||||
!
|
|
||||||
DO ik = 1, nkstot
|
|
||||||
!
|
|
||||||
IF ( lsda ) THEN
|
|
||||||
!
|
|
||||||
IF ( ik == 1 ) WRITE( stdout, 9015 )
|
|
||||||
IF ( ik == ( 1 + nkstot / 2 ) ) WRITE( stdout, 9016 )
|
|
||||||
!
|
|
||||||
END IF
|
|
||||||
!
|
|
||||||
WRITE( stdout, 9020 ) ( xk(i,ik), i = 1, 3 )
|
|
||||||
WRITE( stdout, 9030 ) ( et(ibnd,ik) * rytoev, ibnd = 1, nbnd )
|
|
||||||
!
|
|
||||||
END DO
|
|
||||||
!
|
|
||||||
IF ( lgauss ) THEN
|
|
||||||
!
|
|
||||||
call efermig (et, nbnd, nks, nelec, wk, degauss, ngauss, ef, 0, isk)
|
|
||||||
WRITE( stdout, 9040 ) ef * rytoev
|
|
||||||
!
|
|
||||||
ELSE IF ( ltetra ) THEN
|
|
||||||
!
|
|
||||||
CALL efermit (et, nbnd, nks, nelec, nspin, ntetra, tetra, ef, 0, isk)
|
|
||||||
WRITE( stdout, 9040 ) ef * rytoev
|
|
||||||
!
|
|
||||||
END IF
|
|
||||||
!
|
|
||||||
#if defined (FLUSH)
|
|
||||||
CALL flush( stdout )
|
|
||||||
#endif
|
|
||||||
!
|
|
||||||
! ... do a Berry phase polarization calculation if required
|
|
||||||
!
|
|
||||||
IF ( lberry ) CALL c_phase()
|
|
||||||
!
|
|
||||||
IF ( output_drho /= ' ' ) CALL remove_atomic_rho()
|
|
||||||
!
|
|
||||||
CALL stop_clock( 'electrons' )
|
|
||||||
!
|
|
||||||
RETURN
|
RETURN
|
||||||
!
|
!
|
||||||
END IF
|
END IF
|
||||||
|
@ -296,6 +232,18 @@ SUBROUTINE electrons()
|
||||||
!
|
!
|
||||||
END IF
|
END IF
|
||||||
!
|
!
|
||||||
|
first = ( iter == 1 )
|
||||||
|
!
|
||||||
|
! ... tr2_min is set to an estimate of the error on the energy
|
||||||
|
! ... due to diagonalization - used only in the first scf iteration
|
||||||
|
!
|
||||||
|
10 CONTINUE
|
||||||
|
IF ( first ) THEN
|
||||||
|
tr2_min = nelec*ethr
|
||||||
|
ELSE
|
||||||
|
tr2_min = 0.D0
|
||||||
|
END IF
|
||||||
|
!
|
||||||
! ... diagonalization of the KS hamiltonian
|
! ... diagonalization of the KS hamiltonian
|
||||||
!
|
!
|
||||||
CALL c_bands( iter, ik_, dr2 )
|
CALL c_bands( iter, ik_, dr2 )
|
||||||
|
@ -309,7 +257,7 @@ SUBROUTINE electrons()
|
||||||
!
|
!
|
||||||
IF ( lda_plus_u ) CALL write_ns()
|
IF ( lda_plus_u ) CALL write_ns()
|
||||||
!
|
!
|
||||||
IF ( iter == 1 .AND. lda_plus_u .AND. &
|
IF ( first .AND. lda_plus_u .AND. &
|
||||||
startingpot == 'atomic' .AND. istep == 1 ) CALL ns_adj()
|
startingpot == 'atomic' .AND. istep == 1 ) CALL ns_adj()
|
||||||
!
|
!
|
||||||
! ... calculate total and absolute magnetization
|
! ... calculate total and absolute magnetization
|
||||||
|
@ -331,22 +279,6 @@ SUBROUTINE electrons()
|
||||||
!
|
!
|
||||||
END IF
|
END IF
|
||||||
!
|
!
|
||||||
IF ( iter == 1 ) THEN
|
|
||||||
!
|
|
||||||
! ... for the first scf iteration ethr_min is set for a check
|
|
||||||
! ... in mix_rho ( in mix_rho ethr_min = dr2 * ethr_min )
|
|
||||||
!
|
|
||||||
ethr_min = 1.D0 / nelec
|
|
||||||
!
|
|
||||||
ELSE
|
|
||||||
!
|
|
||||||
! ... otherwise ethr_min is set to a negative number:
|
|
||||||
! ... no check is needed
|
|
||||||
!
|
|
||||||
ethr_min = - 1.D0
|
|
||||||
!
|
|
||||||
END IF
|
|
||||||
!
|
|
||||||
IF ( noncolin ) THEN
|
IF ( noncolin ) THEN
|
||||||
!
|
!
|
||||||
CALL mix_rho_nc( rho, rho_save, nsnew, ns, mixing_beta, dr2, iter, &
|
CALL mix_rho_nc( rho, rho_save, nsnew, ns, mixing_beta, dr2, iter, &
|
||||||
|
@ -354,61 +286,28 @@ SUBROUTINE electrons()
|
||||||
!
|
!
|
||||||
ELSE
|
ELSE
|
||||||
!
|
!
|
||||||
CALL mix_rho( rho, rho_save, nsnew, ns, mixing_beta, dr2, ethr, &
|
CALL mix_rho( rho, rho_save, nsnew, ns, mixing_beta, dr2, &
|
||||||
ethr_min, iter, nmix, flmix, conv_elec )
|
tr2_min, iter, nmix, flmix, conv_elec )
|
||||||
!
|
!
|
||||||
END IF
|
END IF
|
||||||
!
|
!
|
||||||
! ... for the first scf iteration it is controlled that the threshold
|
! ... for the first scf iteration it is controlled that the threshold
|
||||||
! ... is small enough for the diagonalization to be adequate
|
! ... is small enough for the diagonalization to be adequate
|
||||||
!
|
!
|
||||||
IF ( iter == 1 .AND. ethr >= ethr_min ) THEN
|
IF ( first ) THEN
|
||||||
!
|
!
|
||||||
! ... a new diagonalization is needed
|
first = .FALSE.
|
||||||
!
|
IF ( dr2 < tr2_min ) THEN
|
||||||
WRITE( stdout, '(/,5X,"Threshold (ethr) on eigenvalues was ", &
|
|
||||||
& "too large:",/, &
|
|
||||||
& 5X,"Diagonalizing with lowered threshold",/)' )
|
|
||||||
!
|
|
||||||
CALL c_bands( iter, ik_, dr2 )
|
|
||||||
!
|
|
||||||
! ... check if the maximum CPU time has been exceeded
|
|
||||||
! ... or if the user has required a soft exit
|
|
||||||
!
|
|
||||||
IF ( check_stop_now() ) RETURN
|
|
||||||
!
|
|
||||||
CALL sum_band()
|
|
||||||
!
|
|
||||||
IF ( lda_plus_u ) CALL write_ns()
|
|
||||||
!
|
|
||||||
! ... calculate total and absolute magnetization
|
|
||||||
!
|
|
||||||
IF ( lsda .OR. noncolin ) CALL compute_magnetization()
|
|
||||||
!
|
|
||||||
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, vnew )
|
|
||||||
!
|
|
||||||
deband = delta_e ( )
|
|
||||||
!
|
|
||||||
IF ( lda_plus_u .AND. iter <= niter_with_fixed_ns ) THEN
|
|
||||||
!
|
!
|
||||||
ldim2 = ( 2 * Hubbard_lmax + 1 )**2
|
! ... a new diagonalization is needed
|
||||||
nsnew = ns
|
|
||||||
!
|
!
|
||||||
END IF
|
WRITE( stdout, '(/,5X,"Threshold (ethr) on eigenvalues was ", &
|
||||||
!
|
& "too large:",/, &
|
||||||
! ... ethr_min is set to a negative number: no check is needed
|
& 5X,"Diagonalizing with lowered threshold",/)' )
|
||||||
!
|
|
||||||
IF ( noncolin ) THEN
|
|
||||||
!
|
!
|
||||||
CALL mix_rho_nc( rho, rho_save, nsnew, ns, mixing_beta, dr2, &
|
ethr = dr2 / nelec
|
||||||
iter, nmix, flmix, conv_elec )
|
|
||||||
!
|
!
|
||||||
ELSE
|
GO TO 10
|
||||||
!
|
|
||||||
CALL mix_rho( rho, rho_save, nsnew, ns, mixing_beta, dr2, ethr, &
|
|
||||||
-1.D0, iter, nmix, flmix, conv_elec )
|
|
||||||
!
|
!
|
||||||
END IF
|
END IF
|
||||||
!
|
!
|
||||||
|
@ -719,7 +618,6 @@ SUBROUTINE electrons()
|
||||||
!
|
!
|
||||||
9000 FORMAT(/' total cpu time spent up to now is ',F9.2,' secs')
|
9000 FORMAT(/' total cpu time spent up to now is ',F9.2,' secs')
|
||||||
9001 FORMAT(/' Self-consistent Calculation')
|
9001 FORMAT(/' Self-consistent Calculation')
|
||||||
9002 FORMAT(/' Band Structure Calculation')
|
|
||||||
9010 FORMAT(/' iteration #',I3,' ecut=',F9.2,' ryd',5X,'beta=',F4.2)
|
9010 FORMAT(/' iteration #',I3,' ecut=',F9.2,' ryd',5X,'beta=',F4.2)
|
||||||
9015 FORMAT(/' ------ SPIN UP ------------'/)
|
9015 FORMAT(/' ------ SPIN UP ------------'/)
|
||||||
9016 FORMAT(/' ------ SPIN DOWN ----------'/)
|
9016 FORMAT(/' ------ SPIN DOWN ----------'/)
|
||||||
|
@ -759,14 +657,94 @@ SUBROUTINE electrons()
|
||||||
/' potential mean squ. error =',1PE15.1,' ryd^2')
|
/' potential mean squ. error =',1PE15.1,' ryd^2')
|
||||||
9086 FORMAT(/'! total energy =',0PF15.8,' ryd' &
|
9086 FORMAT(/'! total energy =',0PF15.8,' ryd' &
|
||||||
/' potential mean squ. error =',1PE15.1,' ryd^2')
|
/' potential mean squ. error =',1PE15.1,' ryd^2')
|
||||||
! 9090 FORMAT(/' the final potential is written on file ',A14)
|
|
||||||
! 9100 FORMAT(/' this iteration took ',F9.2,' cpu secs')
|
|
||||||
9101 FORMAT(/' End of self-consistent calculation')
|
9101 FORMAT(/' End of self-consistent calculation')
|
||||||
9102 FORMAT(/' End of band structure calculation')
|
|
||||||
9110 FORMAT(/' convergence has been achieved')
|
9110 FORMAT(/' convergence has been achieved')
|
||||||
9120 FORMAT(/' convergence NOT achieved, stopping')
|
9120 FORMAT(/' convergence NOT achieved, stopping')
|
||||||
!
|
!
|
||||||
CONTAINS
|
CONTAINS
|
||||||
|
!
|
||||||
|
!-----------------------------------------------------------------------
|
||||||
|
SUBROUTINE non_scf()
|
||||||
|
!-----------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
USE bp, ONLY : lberry
|
||||||
|
!
|
||||||
|
IMPLICIT NONE
|
||||||
|
!
|
||||||
|
WRITE( stdout, 9002 )
|
||||||
|
!
|
||||||
|
#if defined (FLUSH)
|
||||||
|
CALL flush( stdout )
|
||||||
|
#endif
|
||||||
|
!
|
||||||
|
IF ( imix >= 0 ) rho_save = rho
|
||||||
|
!
|
||||||
|
iter = 1
|
||||||
|
!
|
||||||
|
! ... diagonalization of the KS hamiltonian
|
||||||
|
!
|
||||||
|
CALL c_bands( iter, ik_, dr2 )
|
||||||
|
!
|
||||||
|
conv_elec = .TRUE.
|
||||||
|
!
|
||||||
|
CALL poolrecover( et, nbnd, nkstot, nks )
|
||||||
|
!
|
||||||
|
tcpu = get_clock( 'PWSCF' )
|
||||||
|
WRITE( stdout, 9000 ) tcpu
|
||||||
|
!
|
||||||
|
WRITE( stdout, 9102 )
|
||||||
|
!
|
||||||
|
! ... write band eigenvalues
|
||||||
|
!
|
||||||
|
DO ik = 1, nkstot
|
||||||
|
!
|
||||||
|
IF ( lsda ) THEN
|
||||||
|
!
|
||||||
|
IF ( ik == 1 ) WRITE( stdout, 9015 )
|
||||||
|
IF ( ik == ( 1 + nkstot / 2 ) ) WRITE( stdout, 9016 )
|
||||||
|
!
|
||||||
|
END IF
|
||||||
|
!
|
||||||
|
WRITE( stdout, 9020 ) ( xk(i,ik), i = 1, 3 )
|
||||||
|
WRITE( stdout, 9030 ) ( et(ibnd,ik) * rytoev, ibnd = 1, nbnd )
|
||||||
|
!
|
||||||
|
END DO
|
||||||
|
!
|
||||||
|
IF ( lgauss ) THEN
|
||||||
|
!
|
||||||
|
call efermig (et, nbnd, nks, nelec, wk, degauss, ngauss, ef, 0, isk)
|
||||||
|
|
||||||
|
WRITE( stdout, 9040 ) ef * rytoev
|
||||||
|
!
|
||||||
|
ELSE IF ( ltetra ) THEN
|
||||||
|
!
|
||||||
|
CALL efermit (et, nbnd, nks, nelec, nspin, ntetra, tetra, ef, 0, isk)
|
||||||
|
WRITE( stdout, 9040 ) ef * rytoev
|
||||||
|
!
|
||||||
|
END IF
|
||||||
|
!
|
||||||
|
#if defined (FLUSH)
|
||||||
|
CALL flush( stdout )
|
||||||
|
#endif
|
||||||
|
!
|
||||||
|
! ... do a Berry phase polarization calculation if required
|
||||||
|
!
|
||||||
|
IF ( lberry ) CALL c_phase()
|
||||||
|
!
|
||||||
|
IF ( output_drho /= ' ' ) CALL remove_atomic_rho()
|
||||||
|
!
|
||||||
|
CALL stop_clock( 'electrons' )
|
||||||
|
!
|
||||||
|
9000 FORMAT(/' total cpu time spent up to now is ',F9.2,' secs')
|
||||||
|
9002 FORMAT(/' Band Structure Calculation')
|
||||||
|
9015 FORMAT(/' ------ SPIN UP ------------'/)
|
||||||
|
9016 FORMAT(/' ------ SPIN DOWN ----------'/)
|
||||||
|
9020 FORMAT(/' k =',3F7.4,' band energies (ev):'/)
|
||||||
|
9030 FORMAT( ' ',8F9.4)
|
||||||
|
9040 FORMAT(/' the Fermi energy is ',F10.4,' ev')
|
||||||
|
9102 FORMAT(/' End of band structure calculation')
|
||||||
|
!
|
||||||
|
END SUBROUTINE non_scf
|
||||||
!
|
!
|
||||||
!-----------------------------------------------------------------------
|
!-----------------------------------------------------------------------
|
||||||
SUBROUTINE compute_magnetization()
|
SUBROUTINE compute_magnetization()
|
||||||
|
|
|
@ -11,7 +11,7 @@
|
||||||
!
|
!
|
||||||
!----------------------------------------------------------------------------
|
!----------------------------------------------------------------------------
|
||||||
SUBROUTINE mix_rho( rhout, rhoin, nsout, nsin, alphamix, &
|
SUBROUTINE mix_rho( rhout, rhoin, nsout, nsin, alphamix, &
|
||||||
dr2, ethr, ethr_min, iter, n_iter, filename, conv )
|
dr2, tr2_min, iter, n_iter, filename, conv )
|
||||||
!----------------------------------------------------------------------------
|
!----------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
! ... Modified Broyden's method for charge density mixing
|
! ... Modified Broyden's method for charge density mixing
|
||||||
|
@ -29,7 +29,6 @@ SUBROUTINE mix_rho( rhout, rhoin, nsout, nsin, alphamix, &
|
||||||
USE control_flags, ONLY : imix, ngm0, tr2
|
USE control_flags, ONLY : imix, ngm0, tr2
|
||||||
USE wvfct, ONLY : gamma_only
|
USE wvfct, ONLY : gamma_only
|
||||||
USE wavefunctions_module, ONLY : psic
|
USE wavefunctions_module, ONLY : psic
|
||||||
USE klist, ONLY : nelec
|
|
||||||
USE parser, ONLY : find_free_unit
|
USE parser, ONLY : find_free_unit
|
||||||
USE cell_base, ONLY : omega
|
USE cell_base, ONLY : omega
|
||||||
!
|
!
|
||||||
|
@ -51,10 +50,8 @@ SUBROUTINE mix_rho( rhout, rhoin, nsout, nsin, alphamix, &
|
||||||
alphamix, &! (in) mixing factor
|
alphamix, &! (in) mixing factor
|
||||||
dr2 ! (out) the estimated errr on the energy
|
dr2 ! (out) the estimated errr on the energy
|
||||||
REAL (KIND=DP) :: &
|
REAL (KIND=DP) :: &
|
||||||
ethr, &! actual threshold for diagonalization
|
tr2_min ! estimated error from diagonalization. If the estimated scf error
|
||||||
ethr_min ! minimal threshold for diagonalization
|
! is smaller than this, exit: a more accurate diagonalization is needed
|
||||||
! if in output ethr >= ethr_min a more accurate
|
|
||||||
! diagonalization is needed
|
|
||||||
LOGICAL :: &
|
LOGICAL :: &
|
||||||
conv ! (out) if true the convergence has been reached
|
conv ! (out) if true the convergence has been reached
|
||||||
INTEGER, PARAMETER :: &
|
INTEGER, PARAMETER :: &
|
||||||
|
@ -141,26 +138,15 @@ SUBROUTINE mix_rho( rhout, rhoin, nsout, nsin, alphamix, &
|
||||||
!
|
!
|
||||||
dr2 = rho_dot_product( rhocout, rhocout ) + ns_dot_product( nsout, nsout )
|
dr2 = rho_dot_product( rhocout, rhocout ) + ns_dot_product( nsout, nsout )
|
||||||
!
|
!
|
||||||
! ... ethr_min is dr2 * 1.D0 / nelec or - 1.0 if no check is required
|
! ... if the self-consistency error (dr2) is smaller than the estimated error due to
|
||||||
|
! diagonalization (tr2_min), exit and leave rhoin and rhout unchanged
|
||||||
!
|
!
|
||||||
IF ( ethr_min > 0.D0 ) THEN
|
IF ( dr2 < tr2_min ) THEN
|
||||||
!
|
!
|
||||||
ethr_min = dr2 * ethr_min
|
DEALLOCATE( rhocin, rhocout )
|
||||||
|
CALL stop_clock( 'mix_rho' )
|
||||||
|
RETURN
|
||||||
!
|
!
|
||||||
IF ( ethr > ethr_min ) THEN
|
|
||||||
!
|
|
||||||
ethr = ethr_min
|
|
||||||
!
|
|
||||||
! ... rhoin and rhout are unchanged
|
|
||||||
!
|
|
||||||
DEALLOCATE( rhocin, rhocout )
|
|
||||||
!
|
|
||||||
CALL stop_clock( 'mix_rho' )
|
|
||||||
!
|
|
||||||
RETURN
|
|
||||||
!
|
|
||||||
END IF
|
|
||||||
!
|
|
||||||
END IF
|
END IF
|
||||||
!
|
!
|
||||||
conv = ( dr2 < tr2 )
|
conv = ( dr2 < tr2 )
|
||||||
|
|
Loading…
Reference in New Issue