From 818ee4b111adbd716bf37793c8140d0d2c9aece0 Mon Sep 17 00:00:00 2001 From: giannozz Date: Sun, 22 May 2005 19:05:09 +0000 Subject: [PATCH] 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 --- PW/electrons.f90 | 242 +++++++++++++++++++++-------------------------- PW/mix_rho.f90 | 32 ++----- 2 files changed, 119 insertions(+), 155 deletions(-) diff --git a/PW/electrons.f90 b/PW/electrons.f90 index 644588ac8..0c15e77a5 100644 --- a/PW/electrons.f90 +++ b/PW/electrons.f90 @@ -57,7 +57,6 @@ SUBROUTINE electrons() niter_with_fixed_ns, Hubbard_lmax, & lda_plus_u USE extfield, ONLY : tefield, etotefield - USE bp, ONLY : lberry USE wavefunctions_module, ONLY : evc, evc_nc USE noncollin_module, ONLY : noncolin, npol, magtot_nc USE noncollin_module, ONLY : factlist, pointlist, pointnum, mcons,& @@ -101,13 +100,12 @@ SUBROUTINE electrons() etotefield_new, &! charge_new ! REAL (KIND=DP) :: & - ethr_min ! minimal threshold for diagonalization at the first scf - ! iteration + tr2_min ! estimated error on the energy coming from diagonalization REAL (KIND=DP), ALLOCATABLE :: & wg_g(:,:) ! temporary array used to recover from pools array wg, ! and then print occupations on stdout LOGICAL :: & - exst + exst, first REAL (KIND=DP) :: dr2old, lambda0 ! ! ... external functions @@ -152,69 +150,7 @@ SUBROUTINE electrons() ! IF ( .NOT. lscf ) THEN ! - WRITE( stdout, 9002 ) - ! -#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' ) - ! + CALL non_scf() RETURN ! END IF @@ -296,6 +232,18 @@ SUBROUTINE electrons() ! 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 ! CALL c_bands( iter, ik_, dr2 ) @@ -309,7 +257,7 @@ SUBROUTINE electrons() ! 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() ! ! ... calculate total and absolute magnetization @@ -331,22 +279,6 @@ SUBROUTINE electrons() ! 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 ! CALL mix_rho_nc( rho, rho_save, nsnew, ns, mixing_beta, dr2, iter, & @@ -354,61 +286,28 @@ SUBROUTINE electrons() ! ELSE ! - CALL mix_rho( rho, rho_save, nsnew, ns, mixing_beta, dr2, ethr, & - ethr_min, iter, nmix, flmix, conv_elec ) + CALL mix_rho( rho, rho_save, nsnew, ns, mixing_beta, dr2, & + tr2_min, iter, nmix, flmix, conv_elec ) ! END IF ! ! ... for the first scf iteration it is controlled that the threshold ! ... 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 - ! - 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 + first = .FALSE. + IF ( dr2 < tr2_min ) THEN ! - ldim2 = ( 2 * Hubbard_lmax + 1 )**2 - nsnew = ns + ! ... a new diagonalization is needed ! - END IF - ! - ! ... ethr_min is set to a negative number: no check is needed - ! - IF ( noncolin ) THEN + WRITE( stdout, '(/,5X,"Threshold (ethr) on eigenvalues was ", & + & "too large:",/, & + & 5X,"Diagonalizing with lowered threshold",/)' ) ! - CALL mix_rho_nc( rho, rho_save, nsnew, ns, mixing_beta, dr2, & - iter, nmix, flmix, conv_elec ) + ethr = dr2 / nelec ! - ELSE - ! - CALL mix_rho( rho, rho_save, nsnew, ns, mixing_beta, dr2, ethr, & - -1.D0, iter, nmix, flmix, conv_elec ) + GO TO 10 ! END IF ! @@ -719,7 +618,6 @@ SUBROUTINE electrons() ! 9000 FORMAT(/' total cpu time spent up to now is ',F9.2,' secs') 9001 FORMAT(/' Self-consistent Calculation') -9002 FORMAT(/' Band Structure Calculation') 9010 FORMAT(/' iteration #',I3,' ecut=',F9.2,' ryd',5X,'beta=',F4.2) 9015 FORMAT(/' ------ SPIN UP ------------'/) 9016 FORMAT(/' ------ SPIN DOWN ----------'/) @@ -759,14 +657,94 @@ SUBROUTINE electrons() /' potential mean squ. error =',1PE15.1,' ryd^2') 9086 FORMAT(/'! total energy =',0PF15.8,' ryd' & /' 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') -9102 FORMAT(/' End of band structure calculation') 9110 FORMAT(/' convergence has been achieved') 9120 FORMAT(/' convergence NOT achieved, stopping') ! 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() diff --git a/PW/mix_rho.f90 b/PW/mix_rho.f90 index a525a5643..f1994d7f1 100644 --- a/PW/mix_rho.f90 +++ b/PW/mix_rho.f90 @@ -11,7 +11,7 @@ ! !---------------------------------------------------------------------------- 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 @@ -29,7 +29,6 @@ SUBROUTINE mix_rho( rhout, rhoin, nsout, nsin, alphamix, & USE control_flags, ONLY : imix, ngm0, tr2 USE wvfct, ONLY : gamma_only USE wavefunctions_module, ONLY : psic - USE klist, ONLY : nelec USE parser, ONLY : find_free_unit USE cell_base, ONLY : omega ! @@ -51,10 +50,8 @@ SUBROUTINE mix_rho( rhout, rhoin, nsout, nsin, alphamix, & alphamix, &! (in) mixing factor dr2 ! (out) the estimated errr on the energy REAL (KIND=DP) :: & - ethr, &! actual threshold for diagonalization - ethr_min ! minimal threshold for diagonalization - ! if in output ethr >= ethr_min a more accurate - ! diagonalization is needed + tr2_min ! estimated error from diagonalization. If the estimated scf error + ! is smaller than this, exit: a more accurate diagonalization is needed LOGICAL :: & conv ! (out) if true the convergence has been reached 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 ) ! - ! ... 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 ! conv = ( dr2 < tr2 )