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

Mixing scheme modified in order to conserve the total charge. 

The charge is renormalized at nelec after it has been generated from atomic
densities. It is also renormalized to nelec after the extrapolation procedure.
Added a further check in electrons. rho2zeta modified so that the forward and
backward procedure are exactly one the opposite of the other. Fixed a bug in
update_pot : the wavefinctions were not extrapolated in parallel runs.
Global cleanup of these routines: level 1 blas removed. C.S.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1519 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
sbraccia 2004-12-21 15:28:01 +00:00
parent 9de8453a2a
commit f22df81366
6 changed files with 474 additions and 484 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
Modules/parser.f90
View File

@ -74,6 +74,37 @@ MODULE parser
!
!
!--------------------------------------------------------------------------
FUNCTION find_free_unit()
!--------------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: find_free_unit
INTEGER :: iunit
LOGICAL :: opnd
!
!
unit_loop: DO iunit = 99, 1, -1
!
INQUIRE( UNIT = iunit, OPENED = opnd )
!
IF ( .NOT. opnd ) THEN
!
find_free_unit = iunit
!
RETURN
!
END IF
!
END DO unit_loop
!
CALL errore( 'find_free_unit()', 'free unit not found ?!?', 1 )
!
RETURN
!
END FUNCTION find_free_unit
!
!--------------------------------------------------------------------------
SUBROUTINE delete_if_present( filename, in_warning )
!--------------------------------------------------------------------------
!
@ -81,37 +112,25 @@ MODULE parser
!
CHARACTER(LEN=*), INTENT(IN) :: filename
LOGICAL, OPTIONAL, INTENT(IN) :: in_warning
LOGICAL :: exst, opnd, warning
LOGICAL :: exst, warning
INTEGER :: iunit
!
INQUIRE( FILE = filename, EXIST = exst )
!
IF ( exst ) THEN
!
unit_loop: DO iunit = 99, 1, - 1
!
INQUIRE( UNIT = iunit, OPENED = opnd )
!
IF ( .NOT. opnd ) THEN
!
warning = .FALSE.
!
IF ( PRESENT( in_warning ) ) warning = in_warning
!
OPEN( UNIT = iunit, FILE = filename , STATUS = 'OLD' )
CLOSE( UNIT = iunit, STATUS = 'DELETE' )
!
IF ( warning ) &
WRITE( UNIT = stdout, FMT = '(/,5X,"WARNING: ",A, &
& " file was present; old file deleted")' ) filename
!
RETURN
!
END IF
!
END DO unit_loop
iunit = find_free_unit()
!
CALL errore( 'delete_if_present', 'free unit not found ?!?', 1 )
warning = .FALSE.
!
IF ( PRESENT( in_warning ) ) warning = in_warning
!
OPEN( UNIT = iunit, FILE = filename , STATUS = 'OLD' )
CLOSE( UNIT = iunit, STATUS = 'DELETE' )
!
IF ( warning ) &
WRITE( UNIT = stdout, FMT = '(/,5X,"WARNING: ",A, &
& " file was present; old file deleted")' ) filename
!
END IF
!

21
PW/electrons.f90
View File

@ -29,7 +29,8 @@ SUBROUTINE electrons()
! ... the separate contributions.
!
USE kinds, ONLY : DP
USE parameters, ONLY : npk
USE parameters, ONLY : npk
USE constants, ONLY : eps8, rytoev
USE io_global, ONLY : stdout, ionode
USE cell_base, ONLY : at, bg, alat, omega, tpiba2
USE ions_base, ONLY : zv, nat, ntyp => nsp, ityp, tau
@ -170,7 +171,7 @@ SUBROUTINE electrons()
END IF
!
WRITE( stdout, 9020 ) ( xk(i,ik), i = 1, 3 )
WRITE( stdout, 9030 ) ( et(ibnd,ik) * 13.6058, ibnd = 1, nbnd )
WRITE( stdout, 9030 ) ( et(ibnd,ik) * rytoev, ibnd = 1, nbnd )
!
END DO
!
@ -475,17 +476,20 @@ SUBROUTINE electrons()
WRITE( stdout, 9020 ) ( xk(i,ik), i = 1, 3 )
END IF
!
WRITE( stdout, 9030 ) ( et(ibnd,ik) * 13.6058, ibnd = 1, nbnd )
WRITE( stdout, 9030 ) ( et(ibnd,ik) * rytoev, ibnd = 1, nbnd )
!
END DO
!
IF ( lgauss .OR. ltetra ) WRITE( stdout, 9040 ) ef * 13.6058
IF ( lgauss .OR. ltetra ) WRITE( stdout, 9040 ) ef * rytoev
!
END IF
!
!
IF ( ( ABS( charge - nelec ) / charge ) > 1.D-7 ) &
WRITE( stdout, 9050 ) charge
!
IF ( ( ABS( charge_new - nelec ) / charge_new ) > 1.D-7 ) &
WRITE( stdout, 9051 ) charge_new
!
etot = eband + ( etxc - etxcc ) + ewld + ehart + deband + demet
!
IF ( lda_plus_u ) etot = etot + eth
@ -497,7 +501,7 @@ SUBROUTINE electrons()
!
IF ( imix >= 0 ) THEN
!
IF ( dr2 > 1.D-8 ) THEN
IF ( dr2 > eps8 ) THEN
!
WRITE( stdout, 9081 ) etot, dr2
!
@ -524,7 +528,7 @@ SUBROUTINE electrons()
!
IF ( imix >= 0 ) THEN
!
IF ( dr2 > 1.D-8 ) THEN
IF ( dr2 > eps8 ) THEN
WRITE( stdout, 9081 ) etot, dr2
ELSE
WRITE( stdout, 9083 ) etot, dr2
@ -540,7 +544,7 @@ SUBROUTINE electrons()
!
IF ( imix >= 0 ) THEN
!
IF ( dr2 > 1.D-8 ) THEN
IF ( dr2 > eps8 ) THEN
!
WRITE( stdout, 9080 ) etot, dr2
!
@ -632,6 +636,7 @@ SUBROUTINE electrons()
9030 FORMAT( ' ',8F9.4)
9040 FORMAT(/' the Fermi energy is ',F10.4,' ev')
9050 FORMAT(/' integrated charge =',F15.8)
9051 FORMAT(/' integrated charge_new =',F15.8)
9060 FORMAT(/' band energy sum =', F15.8,' ryd' &
/' one-electron contribution =', F15.8,' ryd' &
/' hartree contribution =', F15.8,' ryd' &

565
PW/mix_rho.f90
View File

File diff suppressed because it is too large Load Diff

68
PW/potinit.f90
View File

@ -80,7 +80,7 @@ SUBROUTINE potinit()
!
CALL io_pot( -1, TRIM( prefix )//'.rho', rho, nspin )
!
WRITE( stdout, '(/5X,"The initial density is read from file ", A14)' ) &
WRITE( stdout, '(/5X,"The initial density is read from file ",A20)' ) &
TRIM( prefix ) // '.rho'
!
! ... here we compute the potential which correspond to the
@ -89,21 +89,31 @@ SUBROUTINE potinit()
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.0D-4 ) &
WRITE( stdout, '(/5X,"starting charge =",F10.5)') charge
!
!
IF ( ABS( charge - nelec ) / charge > 1.D-7 ) THEN
!
WRITE( stdout, &
'(/,5X,"starting charge =",F10.5)') charge
!
WRITE( stdout, &
'(/,5X,"starting charge renormalised to ",F10.5,/)') nelec
!
rho = rho / charge * nelec
!
END IF
!
ELSE
!
CALL io_pot( -1, TRIM( prefix )//'.pot', vr, nspin )
!
WRITE( stdout, '(/5X,"The initial potential is read from file ", A14)' ) &
WRITE( stdout, &
'(/5X,"The initial potential is read from file ",A20)' ) &
TRIM( prefix ) // '.pot'
!
END IF
!
! ... The occupations ns also need to be read in order to build up
! ... the poten
! ... the potential
!
IF ( lda_plus_u ) THEN
!
@ -121,10 +131,10 @@ SUBROUTINE potinit()
!
END IF
!
CALL reduce( ( ldim * ldim * nspin * nat ), ns )
CALL reduce( ( ldim * ldim * nspin * nat ), ns )
CALL poolreduce( ( ldim * ldim * nspin * nat ), ns )
!
CALL DCOPY( ( ldim * ldim * nspin * nat ), ns, 1, nsnew, 1 )
nsnew = ns
!
END IF
!
@ -145,8 +155,10 @@ SUBROUTINE potinit()
IF ( lda_plus_u ) THEN
!
ldim = 2 * Hubbard_lmax + 1
CALL init_ns
CALL DCOPY( ( ldim * ldim * nspin * nat ), ns, 1, nsnew, 1 )
!
CALL init_ns()
!
nsnew = ns
!
END IF
!
@ -158,11 +170,11 @@ SUBROUTINE potinit()
!
CALL io_pot( -1, input_drho, vr, nspin )
!
WRITE( UNIT = stdout, &
FMT = '(/5X,"a scf correction to at. rho is read from", A14)' ) &
WRITE( stdout, &
'(/5X,"a scf correction to at. rho is read from",A20)' ) &
input_drho
!
CALL DAXPY( nrxx, 1.D0, vr, 1, rho, 1 )
rho = rho + vr
!
END IF
!
@ -172,9 +184,18 @@ SUBROUTINE potinit()
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.0D-4 ) &
WRITE( stdout, '(/5X,"starting charge =",F10.5)') charge
!
IF ( ABS( charge - nelec ) / charge > 1.D-7 ) THEN
!
WRITE( stdout, &
'(/,5X,"starting charge =",F10.5)') charge
!
WRITE( stdout, &
'(/,5X,"starting charge renormalised to ",F10.5,/)') nelec
!
rho = rho / charge * nelec
!
END IF
!
END IF
!
@ -186,18 +207,15 @@ SUBROUTINE potinit()
!
IF ( lda_plus_u ) THEN
!
WRITE( stdout, '(/5X,"Parameters of the lda+U calculation:")')
WRITE( stdout, '(5X,"Number of iteration with fixed ns =",I3)') &
WRITE( stdout, '(/5X,"Parameters of the lda+U calculation:")' )
WRITE( stdout, '(5X,"Number of iteration with fixed ns = ",I3)' ) &
niter_with_fixed_ns
WRITE( stdout, '(5X,"Starting ns and Hubbard U :")')
CALL write_ns
WRITE( stdout, '(5X,"Starting ns and Hubbard U :")' )
!
CALL write_ns()
!
END IF
!
IF ( imix >= 0 ) CALL io_pot( +1, TRIM( prefix )//'.rho', rho, nspin )
!
CALL io_pot( +1, TRIM( prefix )//'.pot', vr, nspin )
!
RETURN
!
END SUBROUTINE potinit

145
PW/rho2zeta.f90
View File

@ -1,71 +1,88 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2004 PWSCF 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,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------------------
subroutine rho2zeta (rho, rho_core, nrxx, nspin, iop)
!--------------------------------------------------------------------
! if (iopi.eq.1) transform the spin up spin down charge density rho(*,is
! into rho(*,1) = ( rho_up + rho_dw ) and
! rho(*,2) = ( rho_up - rho_dw ) / rho_tot = zeta
! if (iopi.eq.-1) do the opposit transformation
!----------------------------------------------------------------------------
SUBROUTINE rho2zeta( rho, rho_core, nrxx, nspin, iop )
!---------------------------------------------------------------------------
!
USE io_global, ONLY : stdout
USE kinds
implicit none
integer :: iop, nspin, nrxx, ir
! the input option
! the number of spin polarizations
! the fft grid dimension
! the counter for fft grid
real(kind=DP) :: rho (nrxx, nspin), rho_core (nrxx), rho_up, rho_dw, &
zeta, rhox
! the scf charge density
! the core charge density
! auxiliary variable for rho up
! auxiliary variable for rho dw
! auxiliary variable for zeta
! auxiliary variable for total rho
real(kind=DP) :: eps
! a small number
parameter (eps = 1.0d-30)
if (nspin.eq.1) return
if (iop.eq.1) then
do ir = 1, nrxx
rhox = rho (ir, 1) + rho (ir, 2) + rho_core (ir)
if (rhox.gt.eps) then
zeta = (rho (ir, 1) - rho (ir, 2) ) / rhox
else
zeta = 0.d0
endif
rho (ir, 1) = rho (ir, 1) + rho (ir, 2)
rho (ir, 2) = zeta
enddo
elseif (iop.eq. - 1) then
do ir = 1, nrxx
rhox = rho (ir, 1) + rho_core (ir)
if (rhox.gt.eps) then
rho_up = 0.5d0 * (rho (ir, 1) + rho (ir, 2) * rhox)
rho_dw = 0.5d0 * (rho (ir, 1) - rho (ir, 2) * rhox)
rho (ir, 1) = rho_up
rho (ir, 2) = rho_dw
else
rho (ir, 1) = 0.d0
rho (ir, 2) = 0.d0
endif
enddo
else
WRITE( stdout , * ) ' iop =', iop
call errore ('mag2zeta', 'wrong iop', 1)
endif
return
end subroutine rho2zeta
! ... if ( iopi == 1 ) transform the spin up spin down charge density
! ... rho(*,is) into :
!
! ... rho(*,1) = ( rho_up + rho_dw ) and
! ... rho(*,2) = ( rho_up - rho_dw ) / rho_tot = zeta
!
! ... if ( iopi == -1) do the opposit transformation
!
USE constants, ONLY : eps32
USE io_global, ONLY : stdout
USE kinds, ONLY : DP
!
IMPLICIT NONE
!
INTEGER :: iop, nspin, nrxx, ir
! the input option
! the number of spin polarizations
! the fft grid dimension
! the counter for fft grid
REAL(kind=DP) :: rho(nrxx,nspin), rho_core(nrxx), &
rho_up, rho_dw, zeta, rhox
! the scf charge density
! the core charge density
! auxiliary variable for rho up
! auxiliary variable for rho dw
! auxiliary variable for zeta
! auxiliary variable for total rho
!
!
IF ( nspin == 1 ) RETURN
!
IF ( iop == 1 ) THEN
!
DO ir = 1, nrxx
!
rhox = rho(ir,1) + rho(ir,2) + rho_core(ir)
!
IF ( rhox > eps32 ) THEN
!
zeta = ( rho(ir,1) - rho(ir,2) ) / rhox
!
ELSE
!
zeta = 0.D0
!
END IF
!
rho(ir,1) = rho(ir,1) + rho(ir,2)
rho(ir,2) = zeta
!
END DO
!
ELSE IF ( iop == - 1 ) THEN
!
DO ir = 1, nrxx
!
rhox = rho(ir,1) + rho_core(ir)
!
rho_up = 0.5D0 * ( rho(ir,1) + rho(ir,2) * rhox )
rho_dw = 0.5D0 * ( rho(ir,1) - rho(ir,2) * rhox )
!
rho(ir,1) = rho_up
rho(ir,2) = rho_dw
!
END DO
!
ELSE
!
WRITE( stdout , '(5X,"iop = ",I5)' ) iop
!
CALL errore( 'mag2zeta', 'wrong iop', 1 )
!
END IF
!
RETURN
!
END SUBROUTINE rho2zeta

92
PW/update_pot.f90
View File

@ -7,6 +7,9 @@
!
#include "f_defs.h"
!
#define ONE (1.D0,0.D0)
#define ZERO (0.D0,0.D0)
!
!----------------------------------------------------------------------------
SUBROUTINE update_pot()
!----------------------------------------------------------------------------
@ -42,12 +45,13 @@ SUBROUTINE update_pot()
!
!
USE control_flags, ONLY : order, history
USE io_files, ONLY : prefix, tmp_dir
USE io_files, ONLY : prefix, tmp_dir, nd_nmbr
USE io_global, ONLY : ionode, ionode_id
USE mp, ONLY : mp_bcast
USE mp_global, ONLY : intra_image_comm
!
IMPLICIT NONE
!
! ... local variables
!
INTEGER :: rho_order, wfc_order
LOGICAL :: exists
!
@ -75,38 +79,36 @@ SUBROUTINE update_pot()
rho_order = MIN( 2, history )
!
INQUIRE( FILE = TRIM( tmp_dir ) // &
& TRIM( prefix ) // '.oldrho2', EXIST = exists )
& TRIM( prefix ) // '.old2rho', EXIST = exists )
!
IF ( exists ) THEN
!
rho_order = MIN( 3, history )
!
END IF
IF ( exists ) rho_order = MIN( 3, history )
!
END IF
!
wfc_order = MIN( 1, history, order )
!
INQUIRE( FILE = TRIM( tmp_dir ) // &
& TRIM( prefix ) // '.oldwfc', EXIST = exists )
!
IF ( exists ) THEN
!
wfc_order = MIN( 2, history, order )
!
INQUIRE( FILE = TRIM( tmp_dir ) // &
& TRIM( prefix ) // '.oldwfc2', EXIST = exists )
!
IF ( exists ) THEN
!
wfc_order = MIN( 3, history, order )
!
END IF
!
END IF
END IF
!
CALL extrapolate_charge( rho_order )
!
IF ( ionode ) THEN
!
wfc_order = MIN( 1, history, order )
!
INQUIRE( FILE = TRIM( tmp_dir ) // &
& TRIM( prefix ) // '.oldwfc' // nd_nmbr, EXIST = exists )
!
IF ( exists ) THEN
!
wfc_order = MIN( 2, history, order )
!
INQUIRE( FILE = TRIM( tmp_dir ) // &
& TRIM( prefix ) // '.old2wfc' // nd_nmbr , EXIST = exists )
!
IF ( exists ) wfc_order = MIN( 3, history, order )
!
END IF
!
END IF
!
CALL mp_bcast( wfc_order, ionode_id, intra_image_comm )
!
IF ( order >= 2 ) CALL extrapolate_wfcs( wfc_order )
!
CALL stop_clock( 'update_pot' )
@ -134,13 +136,12 @@ SUBROUTINE extrapolate_charge( rho_order )
USE cellmd, ONLY : lmovecell, omega_old
USE vlocal, ONLY : strf
USE io_files, ONLY : prefix
USE klist, ONLY : nelec
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: rho_order
!
! ... local variables
!
REAL(KIND=DP), ALLOCATABLE :: work(:), work1(:)
! work is the difference between charge density and atomic charge
! at time t
@ -191,7 +192,7 @@ SUBROUTINE extrapolate_charge( rho_order )
! ... work -> oldrho2
!
CALL io_pot( + 1, TRIM( prefix )//'.oldrho', rho, 1 )
CALL io_pot( + 1, TRIM( prefix )//'.oldrho2', work, 1 )
CALL io_pot( + 1, TRIM( prefix )//'.old2rho', work, 1 )
!
! ... alpha0 has been calculated in move_ions
!
@ -206,14 +207,14 @@ SUBROUTINE extrapolate_charge( rho_order )
! ... oldrho2 -> work1
! ... oldrho -> work
!
CALL io_pot( - 1, TRIM( prefix )//'.oldrho2', work1, 1 )
CALL io_pot( - 1, TRIM( prefix )//'.old2rho', work1, 1 )
CALL io_pot( - 1, TRIM( prefix )//'.oldrho', work, 1 )
!
! ... rho -> oldrho
! ... work -> oldrho2
!
CALL io_pot( + 1, TRIM( prefix )//'.oldrho', rho, 1 )
CALL io_pot( + 1, TRIM( prefix )//'.oldrho2', work, 1 )
CALL io_pot( + 1, TRIM( prefix )//'.old2rho', work, 1 )
!
! ... alpha0 and beta0 have been calculated in move_ions
!
@ -249,6 +250,18 @@ SUBROUTINE extrapolate_charge( rho_order )
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,"extrapolated charge =",F10.5)') charge
!
WRITE( stdout, &
'(/,5X,"extrapolated charge renormalised to ",F10.5,/)') nelec
!
rho = rho / charge * nelec
!
END IF
!
! ... write potential (and rho) on file
!
IF ( imix >= 0 ) CALL io_pot( + 1, TRIM( prefix )//'.rho', rho, nspin )
@ -270,9 +283,6 @@ SUBROUTINE extrapolate_wfcs( wfc_order )
! ... of the basis of the t-dt and t time steps, according to a recipe
! ... by Mead, Rev. Mod. Phys., vol 64, pag. 51 (1992), eqs. 3.20-3.29
!
#define ONE (1.D0,0.D0)
#define ZERO (0.D0,0.D0)
!
USE io_global, ONLY : stdout
USE kinds, ONLY : DP
USE klist, ONLY : nks
@ -286,8 +296,6 @@ SUBROUTINE extrapolate_wfcs( wfc_order )
!
INTEGER, INTENT(IN) :: wfc_order
!
! ... local variables
!
INTEGER :: j, i, ik, zero_ew, lwork, info
! do-loop variables
! counter on k-points
@ -333,7 +341,7 @@ SUBROUTINE extrapolate_wfcs( wfc_order )
CALL diropn( iunoldwfc, TRIM( prefix ) // '.oldwfc', nwordwfc, exst )
!
IF ( order > 2 ) &
CALL diropn( iunoldwfc2, TRIM( prefix ) // '.oldwfc2', nwordwfc, exst )
CALL diropn( iunoldwfc2, TRIM( prefix ) // '.old2wfc', nwordwfc, exst )
!
ALLOCATE( evcold(npwx,nbnd) )
!
@ -450,7 +458,7 @@ SUBROUTINE extrapolate_wfcs( wfc_order )
! ... case : wfc_order = 3
!
CALL diropn( iunoldwfc, TRIM( prefix ) // '.oldwfc', nwordwfc, exst )
CALL diropn( iunoldwfc2, TRIM( prefix ) // '.oldwfc2', nwordwfc, exst )
CALL diropn( iunoldwfc2, TRIM( prefix ) // '.old2wfc', nwordwfc, exst )
!
ALLOCATE( evcold(npwx,nbnd) )
!