scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

- do not set to zero the negative part of the charge density 

when starting from superpositions of atomic charges: it is useless
- print HOMO/top of the VB and LUMO/bottom of the CB when available
- always start from charge density and recalculate the potential,
  also in a non-scf calculation
- do not overwrite the charge density after a non-scf calculation


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1919 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2005-05-27 12:52:50 +00:00
parent 65922c19c0
commit 97dcfedf72
4 changed files with 117 additions and 121 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
PW/atomic_rho.f90
View File

@ -152,41 +152,37 @@ subroutine atomic_rho (rhoa, nspina)
rhoneg = rhoneg + MIN (0.d0, DREAL (psic (ir)) )
rhoima = rhoima + abs (DIMAG (psic (ir) ) )
enddo
rhoneg = rhoneg / (nr1 * nr2 * nr3)
rhoima = rhoima / (nr1 * nr2 * nr3)
rhoneg = omega * rhoneg / (nr1 * nr2 * nr3)
rhoima = omega * rhoima / (nr1 * nr2 * nr3)
#ifdef __PARA
call reduce (1, rhoneg)
call reduce (1, rhoima)
#endif
IF ( rhoima > 1.0d-4 ) THEN
WRITE( stdout,'(/4x," imaginary charge or magnetization ",&
& f12.6," (component ",i1,") set to zero")') rhoima, is
WRITE( stdout,'(/4x," integrated imaginary charge or magnetization:",&
& f12.6," (component ",i1,") set to zero")') rhoima, is
END IF
IF ( (is == 1) .OR. lsda ) THEN
!
! set starting charge rho or rho up or rho down > 0
!
DO ir = 1, nrxx
rhoa (ir, is) = MAX (0.d0, DREAL (psic (ir)) )
END DO
IF ( (rhoneg < -1.0d-4) ) THEN
IF ( lsda ) THEN
WRITE( stdout,'(/4x," negative starting charge (spin=",i2,&
& ") ",f12.6," set to zero")') 2*is-3, rhoneg
WRITE( stdout,'(/4x," integrated negative starting charge ", &
&"(spin=",i2,"):",f12.6)') 2*is-3, rhoneg
ELSE
WRITE( stdout,'(/4x," negative starting charge ", &
& f12.6," set to zero")') rhoneg
WRITE( stdout,'(/4x," integrated negative starting charge: ", &
& f12.6)') rhoneg
END IF
END IF
ELSE
!
! magnetization for non collinear case (is = 2,3,4)
!
DO ir = 1, nrxx
rhoa (ir, is) = DREAL (psic (ir))
END DO
END IF
!
! set imaginary terms to zero - negative terms are not set to zero
! because it is basically useless to do it in real space: negative
! charge will re-appear when Fourier-transformed back and forth
!
DO ir = 1, nrxx
rhoa (ir, is) = DREAL (psic (ir))
END DO
!
enddo
deallocate (rhocg)

40
PW/electrons.f90
View File

@ -75,6 +75,7 @@ SUBROUTINE electrons()
dr2, &! the norm of the diffence between potential
charge, &! the total charge
mag, &! local magnetization
ehomo, elumo, &! highest occupied and lowest onuccupied levels
tcpu ! cpu time
INTEGER :: &
i, &! counter on polarization
@ -432,7 +433,23 @@ SUBROUTINE electrons()
ELSE
WRITE( stdout, 9040 ) ef * rytoev
END IF
END IF
ELSE
!
IF (nspin == 1) THEN
ibnd = nint (nelec) / 2.d0
ELSE
ibnd = nint (nelec)
END IF
!
IF (nbnd > ibnd ) THEN
!
ehomo = MAXVAL ( et( ibnd , 1:nkstot) )
elumo = MINVAL ( et( ibnd+1, 1:nkstot) )
!
WRITE( stdout, 9042 ) ehomo * rytoev, elumo * rytoev
!
END IF
END IF
!
END IF
!
@ -561,6 +578,7 @@ SUBROUTINE electrons()
9021 FORMAT(/' k =',3F7.4,' (',I6,' PWs) bands (ev):'/ )
9030 FORMAT( ' ',8F9.4 )
9032 FORMAT(/' occupation numbers ' )
9042 FORMAT(/' highest occupied, lowest unoccupied level (ev): ',2F10.4 )
9041 FORMAT(/' the spin up/dw Fermi energies are ',2F10.4,' ev' )
9040 FORMAT(/' the Fermi energy is ',F10.4,' ev' )
9050 FORMAT(/' integrated charge =',F15.8 )
@ -604,7 +622,6 @@ SUBROUTINE electrons()
!
IMPLICIT NONE
!
!
WRITE( stdout, 9002 )
!
CALL flush_unit( stdout )
@ -652,6 +669,24 @@ SUBROUTINE electrons()
!
WRITE( stdout, 9040 ) ef * rytoev
!
ELSE
!
IF (nspin == 1) THEN
ibnd = nint (nelec) / 2.d0
ELSE
ibnd = nint (nelec)
END IF
!
IF (nbnd > ibnd ) THEN
!
ehomo = MAXVAL ( et( ibnd , 1:nkstot) )
elumo = MINVAL ( et( ibnd+1, 1:nkstot) )
!
IF ( ehomo < elumo ) &
WRITE( stdout, 9042 ) ehomo * rytoev, elumo * rytoev
!
END IF
!
END IF
!
CALL flush_unit( stdout )
@ -671,6 +706,7 @@ SUBROUTINE electrons()
9020 FORMAT(/' k =',3F7.4,' band energies (ev):'/ )
9030 FORMAT( ' ',8F9.4 )
9040 FORMAT(/' the Fermi energy is ',F10.4,' ev' )
9042 FORMAT(/' Highest occupied, lowest unoccupied level (ev): ',2F10.4 )
9102 FORMAT(/' End of band structure calculation' )
!
END SUBROUTINE non_scf

133
PW/potinit.f90
View File

@ -14,11 +14,14 @@ SUBROUTINE potinit()
! ... This routine initializes the self consistent potential in the array
! ... vr. There are three possible cases:
!
! ... a) In this run the code is restarting from a broken run
! ... b) The potential (or rho) is read from file
! ... c) if a and b are both false, the total charge is computed
! ... as a sum of atomic charges, and the corresponding potential
! ... is saved in vr
! ... a) the code is restarting from a broken run:
! ... read rho from data stored during the previous run
! ... b) the code is performing a non-scf calculation following a scf one:
! ... read rho from the file produced by the scf calculation
! ... c) the code starts a new calculation:
! ... calculate rho as a sum of atomic charges
!
! ... In all cases the scf potential is recalculated and saved in vr
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
@ -55,15 +58,7 @@ SUBROUTINE potinit()
!
IF ( ionode ) THEN
!
IF ( lscf ) THEN
!
CALL seqopn( 4, TRIM( prefix )//'.rho', 'UNFORMATTED', exst )
!
ELSE
!
CALL seqopn( 4, TRIM( prefix )//'.pot', 'UNFORMATTED', exst )
!
END IF
CALL seqopn( 4, TRIM( prefix )//'.rho', 'UNFORMATTED', exst )
!
IF ( exst ) THEN
!
@ -81,49 +76,16 @@ SUBROUTINE potinit()
!
IF ( startingpot == 'file' .AND. exst ) THEN
!
! ... First case, the potential is read from file
! ... NB: this case applies also for a restarting run, in which case
! ... potential and rho files have been read from the restart file
! ... Cases a) and b): the charge density is read from file
!
CALL io_pot( -1, TRIM( prefix )//'.rho', rho, nspin )
!
IF ( lscf ) THEN
!
CALL io_pot( -1, TRIM( prefix )//'.rho', rho, nspin )
!
WRITE( stdout, '(/5X,"The initial density is read from file ", A20)' ) &
TRIM( prefix ) // '.rho'
!
! ... here we compute the potential which correspond to the
! ... initial charge
!
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 )
!
!
IF ( ABS( charge - nelec ) / charge > 1.D-7 ) THEN
!
WRITE( stdout, &
'(/,5X,"starting charge ",F10.5,", renormalised to ",F10.5)') &
charge, nelec
!
rho = rho / charge * nelec
!
! ... and compute v_of_rho again
!
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 )
!
END IF
!
WRITE( stdout, '(/5X,"The initial density is read from file ", A20)' )&
TRIM( prefix ) // '.rho'
ELSE
!
CALL io_pot( -1, TRIM( prefix )//'.pot', vr, nspin )
!
WRITE( stdout, &
'(/5X,"The initial potential is read from file ", A20)' ) &
TRIM( prefix ) // '.pot'
!
WRITE( stdout, '(/5X,"The potential is recalculated from file ",A20)')&
TRIM( prefix ) // '.rho'
END IF
!
! ... The occupations ns also need to be read in order to build up
@ -154,11 +116,11 @@ SUBROUTINE potinit()
!
ELSE
!
! ... Second case, the potential is built from a superposition
! ... Case c): the potential is built from a superposition
! ... of atomic charges contained in the array rho_at
!
!
IF ( startingpot == 'file' .AND. .NOT. exst ) &
WRITE( stdout, '(5X,"Cannot read pot/rho file: not found")' )
WRITE( stdout, '(5X,"Cannot read rho : file not found")' )
!
WRITE( UNIT = stdout, &
FMT = '(/5X,"Initial potential from superposition of free atoms")' )
@ -166,8 +128,6 @@ SUBROUTINE potinit()
! ... in the lda+U case set the initial value of ns
!
IF ( lda_plus_u ) THEN
!
ldim = 2 * Hubbard_lmax + 1
!
CALL init_ns()
!
@ -191,32 +151,41 @@ SUBROUTINE potinit()
!
END IF
!
! ... here we compute the potential which corresponds to the
! ... initial charge
!
!
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 )
!
IF ( ABS( charge - nelec ) / charge > 1.D-7 ) THEN
!
WRITE( stdout, &
'(/,5X,"starting charge ",F10.5,", renormalised to ",F10.5)') &
charge, nelec
!
rho = rho / charge * nelec
!
! ... and compute v_of_rho again
!
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 )
!
END IF
!
END IF
!
! ... check the integral of the starting charge
!
IF ( nspin == 2 ) THEN
!
charge = SUM ( rho (:, 1:nspin) ) * omega / ( nr1 * nr2 * nr3 )
!
ELSE
!
charge = SUM ( rho (:, 1) ) * omega / ( nr1 * nr2 * nr3 )
!
END IF
!
IF ( lscf .AND. ABS( charge - nelec ) / charge > 1.D-6 ) THEN
!
WRITE( stdout, &
'(/,5X,"starting charge ",F10.5,", renormalised to ",F10.5)') &
charge, nelec
!
rho = rho / charge * nelec
!
ELSE IF ( .NOT. lscf .AND. ABS( charge - nelec ) / charge > 1.D-6 ) THEN
!
CALL errore ( 'potinit', 'starting and expected charges differ', 1 )
!
END IF
!
! ... compute the potential and store it in vr
!
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 )
!
! ... define the total local potential (external+scf)
!
CALL set_vrs( vrs, vltot, vr, nrxx, nspin, doublegrid )

27
PW/punch.f90
View File

@ -10,12 +10,10 @@ SUBROUTINE punch()
!-----------------------------------------------------------------------
!
! ... This routine is called at the end of the run to save on a file
! ... the information needed to the phonon program.
! ... the information needed for further processing
!
USE io_global, ONLY : stdout
USE klist, ONLY : nks, nkstot
USE lsda_mod, ONLY : nspin
USE scf, ONLY : rho
USE control_flags, ONLY : reduce_io, lscf
USE wvfct, ONLY : et, wg, nbnd
USE wavefunctions_module, ONLY : evc, evc_nc
@ -30,8 +28,7 @@ SUBROUTINE punch()
LOGICAL :: exst
!
!
WRITE( UNIT = stdout, &
FMT = '(/,5X,"Writing file ",A14," for program phonon")' ) &
WRITE( UNIT = stdout, FMT = '(/,5X,"Writing data file ",A14)' ) &
TRIM( prefix ) // '.save'
!
kunittmp = 1
@ -50,11 +47,14 @@ SUBROUTINE punch()
!
ENDIF
!
! ... The following instruction is used when more k-points are needed
! ... for finite-q phonon calculations (on fine q-grid) then those needed
! ... for self-consistency. In such a case, a self-consistent calculation
! ... with few k-points is followed by a non-self-consistent one with added
! ... k-points, whose weight is set to zero.
! ... The following instruction is needed to recalculate the weights
! ... of k-points: this is useful for finite-q phonon calculations
! ... and when more k-points are needed than in self-consistency.
! ... In such a case, a self-consistent calculation with few k-points
! ... is followed by a non-self-consistent one with added k-points,
! ... whose weight is set to zero. Note that the charge density
! ... is recalculated but NOT written to file, because doing this
! ... might spoil the charge density in other cases
!
IF ( .NOT. lscf ) CALL sum_band()
!
@ -66,7 +66,7 @@ SUBROUTINE punch()
! ... xk, wk, isk, et, wg are distributed across pools
! ... the first node has a complete copy of xk, wk, isk,
! ... while eigenvalues et and weights wg must be
! ... explicitely collected to the first node
! ... explicitly collected to the first node
!
CALL poolrecover( et, nbnd, nkstot, nks )
CALL poolrecover( wg, nbnd, nkstot, nks )
@ -76,13 +76,8 @@ SUBROUTINE punch()
kunittmp = kunit
!
#endif
!
! ... Write the charge density on a separate file
!
CALL io_pot( + 1, TRIM(prefix)//'.rho', rho, nspin )
!
iunpun = 4
!
CALL writefile_new( 'all', iunpun, et, wg, kunittmp )
!
RETURN