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

Split ef_shift and ef_shift_wfc

This commit is contained in:
Jae-Mo Lihm 2021-09-08 13:41:52 +09:00
parent eed27240e3
commit a3246a258f
6 changed files with 169 additions and 132 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
HP/src/hp_print_clock.f90
View File

@ -37,7 +37,7 @@ SUBROUTINE hp_print_clock
CALL print_clock ('hp_calc_chi')
CALL print_clock ('hp_postproc')
CALL print_clock ('hp_vpsifft')
CALL print_clock ('hp_ef_shift')
CALL print_clock ('ef_shift')
CALL print_clock ('hp_run_nscf')
CALL print_clock ('hp_postproc')
!

4
HP/src/hp_solve_linear_system.f90
View File

@ -260,9 +260,9 @@ SUBROUTINE hp_solve_linear_system (na, iq)
IF (lmetq0) THEN
!
IF (okpaw) THEN
CALL ef_shift(.FALSE., 1, dos_ef, ldos, ldoss, drhoscfh, dbecsum=dbecsum, becsum1=becsum1)
CALL ef_shift(1, dos_ef, ldos, drhoscfh, dbecsum=dbecsum, becsum1=becsum1)
ELSE
CALL ef_shift(.FALSE., 1, dos_ef, ldos, ldoss, drhoscfh)
CALL ef_shift(1, dos_ef, ldos, drhoscfh)
ENDIF
!
! Check that def is not too large (it is in Ry).

2
LR_Modules/Makefile
View File

@ -68,7 +68,7 @@ lanczos_nonhermitian.o \
lanczos_pseudohermitian_c.o \
lanczos_nonhermitian_c.o \
response_kernels.o \
efermi_shift.o \
efermi_shift.o
TLDEPS=mods pwlibs

273
LR_Modules/efermi_shift.f90
View File

@ -11,7 +11,7 @@ MODULE efermi_shift
! the change of the Fermi energy for each pert.
! NB: def(3) should be def (npertx), but it is used only at Gamma
! where the dimension of irreps never exceeds 3
!
! Define an abstract interface to use a callback
ABSTRACT INTERFACE
SUBROUTINE def_symmetrization(def, irr)
@ -20,12 +20,11 @@ MODULE efermi_shift
COMPLEX(DP) :: def(3)
END SUBROUTINE
END INTERFACE
CONTAINS
!
CONTAINS
!
!-----------------------------------------------------------------------
SUBROUTINE ef_shift (update_wfc, npert, dos_ef, ldos, ldoss, drhoscf, &
dbecsum, becsum1, irr, sym_def)
SUBROUTINE ef_shift (npert, dos_ef, ldos, drhoscf, dbecsum, becsum1, irr, sym_def)
!-----------------------------------------------------------------------
!! This routine takes care of the effects of a shift of Ef, due to the
!! perturbation, that can take place in a metal at q=0
@ -36,38 +35,24 @@ SUBROUTINE ef_shift (update_wfc, npert, dos_ef, ldos, ldoss, drhoscf, &
USE mp, ONLY : mp_sum
USE io_global, ONLY : stdout
USE ions_base, ONLY : nat
USE wavefunctions, ONLY : evc
USE cell_base, ONLY : omega
USE fft_base, ONLY : dfftp, dffts
USE fft_base, ONLY : dfftp
USE fft_interfaces, ONLY : fwfft, invfft
USE gvect, ONLY : gg
USE buffers, ONLY : get_buffer, save_buffer
USE lsda_mod, ONLY : nspin
USE uspp_param, ONLY : nhm
USE wvfct, ONLY : npwx, et
USE klist, ONLY : degauss, ngauss, ngk, ltetra
USE ener, ONLY : ef
USE noncollin_module, ONLY : noncolin, npol, nspin_mag, nspin_lsda
USE qpoint, ONLY : nksq
USE control_lr, ONLY : nbnd_occ
USE units_lr, ONLY : iuwfc, lrwfc, lrdwf, iudwf
USE eqv, ONLY : dpsi
USE dfpt_tetra_mod, ONLY : dfpt_tetra_delta
USE noncollin_module, ONLY : nspin_mag, nspin_lsda
!
IMPLICIT NONE
!
! input/output variables
!
LOGICAL, INTENT(IN) :: update_wfc
!! if true the eigenfunctions are updated
INTEGER, INTENT(IN) :: npert
!! the number of perturbation
REAL(DP), INTENT(IN) :: dos_ef
!! density of states at Ef
COMPLEX(DP), INTENT(IN) :: ldos(dfftp%nnr, nspin_mag)
!! local DOS at Ef (with augmentation)
COMPLEX(DP), INTENT(IN) :: ldoss(dffts%nnr, nspin_mag)
!! local DOS at Ef without augmentation
COMPLEX(DP), INTENT(INOUT) :: drhoscf(dfftp%nnr, nspin_mag, npert)
!! the change of the charge (with augmentation)
COMPLEX(DP), INTENT(INOUT), OPTIONAL :: dbecsum((nhm*(nhm+1))/2, nat, nspin_mag, npert)
@ -83,16 +68,101 @@ SUBROUTINE ef_shift (update_wfc, npert, dos_ef, ldos, ldoss, drhoscf, &
!
! local variables
!
!--> these quantities may be complex since perturbation may be
complex(DP) :: delta_n, wfshift
! the change in electron number
! the shift coefficient for the wavefunction
real(DP), external :: w0gauss
! the smeared delta function
integer :: npw, ibnd, ik, is, ipert, nrec, ikrec
INTEGER :: is
!! counter on spin polarizations
INTEGER :: ipert
!! counter on perturbations
COMPLEX(DP) :: delta_n
!! the change in electron number
!! This may be complex since perturbation may be complex
!
REAL(DP), external :: w0gauss
!! the smeared delta function
!
call start_clock ('ef_shift')
!
IF (npert > 3) CALL errore("ef_shift", "npert exceeds 3", 1)
!
! determines Fermi energy shift (such that each pertubation is neutral)
!
WRITE( stdout, * )
do ipert = 1, npert
delta_n = (0.d0, 0.d0)
do is = 1, nspin_lsda
CALL fwfft ('Rho', drhoscf(:,is,ipert), dfftp)
if (gg(1) < 1.0d-8) delta_n = delta_n + omega*drhoscf(dfftp%nl(1),is,ipert)
CALL invfft ('Rho', drhoscf(:,is,ipert), dfftp)
enddo
call mp_sum ( delta_n, intra_bgrp_comm )
IF ( ABS(dos_ef) > 1.d-18 ) THEN
def (ipert) = - delta_n / dos_ef
ELSE
def (ipert) = 0.0_dp
ENDIF
enddo
!
! symmetrizes the Fermi energy shift
!
IF (present(sym_def)) CALL sym_def(def, irr)
WRITE( stdout, '(5x,"Pert. #",i3,": Fermi energy shift (Ry) =",2es15.4)')&
(ipert, def (ipert) , ipert = 1, npert )
!
! corrects the density response accordingly...
!
do ipert = 1, npert
call zaxpy (dfftp%nnr*nspin_mag, def(ipert), ldos, 1, drhoscf(1,1,ipert), 1)
enddo
!
! In the PAW case there is also a metallic term
!
IF (PRESENT(dbecsum) .AND. PRESENT(becsum1)) THEN
DO ipert = 1, npert
dbecsum(:,:,:,ipert) = dbecsum(:,:,:,ipert) &
+ def(ipert) * CMPLX(becsum1(:,:,:)*0.5_DP, 0.0_DP, KIND=DP)
ENDDO
ENDIF
!
CALL stop_clock ('ef_shift')
!
END SUBROUTINE ef_shift
!-------------------------------------------------------------------------
!
!-------------------------------------------------------------------------
SUBROUTINE ef_shift_wfc(npert, ldoss, drhoscf)
!-----------------------------------------------------------------------
!! This routine takes care of the effects of a shift of Ef, due to the
!! perturbation, that can take place in a metal at q=0, on the wavefunction
!
USE kinds, ONLY : DP
USE mp, ONLY : mp_sum
USE wavefunctions, ONLY : evc
USE fft_base, ONLY : dfftp, dffts
USE fft_interfaces, ONLY : fwfft, invfft
USE buffers, ONLY : get_buffer, save_buffer
USE wvfct, ONLY : npwx, et
USE klist, ONLY : degauss, ngauss, ngk, ltetra
USE ener, ONLY : ef
USE noncollin_module, ONLY : noncolin, npol, nspin_mag
USE qpoint, ONLY : nksq
USE control_lr, ONLY : nbnd_occ
USE units_lr, ONLY : iuwfc, lrwfc, lrdwf, iudwf
USE eqv, ONLY : dpsi
USE dfpt_tetra_mod, ONLY : dfpt_tetra_delta
!
IMPLICIT NONE
!
! input/output variables
!
INTEGER, INTENT(IN) :: npert
!! the number of perturbation
COMPLEX(DP), INTENT(IN) :: ldoss(dffts%nnr, nspin_mag)
!! local DOS at Ef without augmentation
COMPLEX(DP), INTENT(INOUT) :: drhoscf(dfftp%nnr, nspin_mag, npert)
!! the change of the charge (with augmentation)
!
! local variables
!
INTEGER :: npw, ibnd, ik, is, ipert, nrec, ikrec
! counter on occupied bands
! counter on k-point
! counter on spin polarizations
@ -100,96 +170,63 @@ SUBROUTINE ef_shift (update_wfc, npert, dos_ef, ldos, ldoss, drhoscf, &
! record number
! record position of wfc at k
! auxiliary for spin
COMPLEX(DP) :: wfshift
!! the shift coefficient for the wavefunction
!! This may be complex since perturbation may be complex
!
call start_clock ('ef_shift')
REAL(DP), external :: w0gauss
! the smeared delta function
!
call start_clock ('ef_shift_wfc')
!
IF (npert > 3) CALL errore("ef_shift", "npert exceeds 3", 1)
!
if (.not.update_wfc) then
! determines Fermi energy shift (such that each pertubation is neutral)
WRITE( stdout, * )
do ipert = 1, npert
delta_n = (0.d0, 0.d0)
do is = 1, nspin_lsda
CALL fwfft ('Rho', drhoscf(:,is,ipert), dfftp)
if (gg(1) < 1.0d-8) delta_n = delta_n + omega*drhoscf(dfftp%nl(1),is,ipert)
CALL invfft ('Rho', drhoscf(:,is,ipert), dfftp)
enddo
call mp_sum ( delta_n, intra_bgrp_comm )
IF ( ABS(dos_ef) > 1.d-18 ) THEN
def (ipert) = - delta_n / dos_ef
ELSE
def (ipert) = 0.0_dp
ENDIF
enddo
!
! symmetrizes the Fermi energy shift
!
IF (present(sym_def)) CALL sym_def(def, irr)
WRITE( stdout, '(5x,"Pert. #",i3,": Fermi energy shift (Ry) =",2es15.4)')&
(ipert, def (ipert) , ipert = 1, npert )
!
! corrects the density response accordingly...
!
do ipert = 1, npert
call zaxpy (dfftp%nnr*nspin_mag, def(ipert), ldos, 1, drhoscf(1,1,ipert), 1)
enddo
!
! In the PAW case there is also a metallic term
!
IF (PRESENT(dbecsum) .AND. PRESENT(becsum1)) THEN
DO ipert = 1, npert
dbecsum(:,:,:,ipert) = dbecsum(:,:,:,ipert) &
+ def(ipert) * CMPLX(becsum1(:,:,:)*0.5_DP, 0.0_DP, KIND=DP)
ENDDO
ENDIF
!
else
!
! does the same for perturbed wfc
!
do ik = 1, nksq
npw = ngk (ik)
!
! reads unperturbed wavefuctions psi_k in G_space, for all bands
!
ikrec = ik
if (nksq > 1) call get_buffer (evc, lrwfc, iuwfc, ikrec)
!
! reads delta_psi from iunit iudwf, k=kpoint
!
do ipert = 1, npert
nrec = (ipert - 1) * nksq + ik
IF (nksq > 1 .OR. npert > 1) CALL get_buffer(dpsi, lrdwf, iudwf, nrec)
do ibnd = 1, nbnd_occ (ik)
!
if(ltetra) then
wfshift = 0.5d0 * def(ipert) * dfpt_tetra_delta(ibnd,ik)
else
wfshift = 0.5d0 * def(ipert) * w0gauss( (ef-et(ibnd,ik))/degauss, ngauss) / degauss
end if
!
IF (noncolin) THEN
call zaxpy (npwx*npol,wfshift,evc(1,ibnd),1,dpsi(1,ibnd),1)
ELSE
call zaxpy (npw, wfshift, evc(1,ibnd), 1, dpsi(1,ibnd), 1)
ENDIF
enddo
!
! writes corrected delta_psi to iunit iudwf, k=kpoint,
!
IF (nksq > 1 .OR. npert > 1) CALL save_buffer(dpsi, lrdwf, iudwf, nrec)
enddo
enddo
do ipert = 1, npert
do is = 1, nspin_mag
call zaxpy (dffts%nnr, def(ipert), ldoss(1,is), 1, drhoscf(1,is,ipert), 1)
enddo
enddo
endif
! Update the perturbed wavefunctions according to the Fermi energy shift
!
CALL stop_clock ('ef_shift')
do ik = 1, nksq
npw = ngk (ik)
!
! reads unperturbed wavefuctions psi_k in G_space, for all bands
!
ikrec = ik
if (nksq > 1) call get_buffer (evc, lrwfc, iuwfc, ikrec)
!
! reads delta_psi from iunit iudwf, k=kpoint
!
do ipert = 1, npert
nrec = (ipert - 1) * nksq + ik
IF (nksq > 1 .OR. npert > 1) CALL get_buffer(dpsi, lrdwf, iudwf, nrec)
do ibnd = 1, nbnd_occ (ik)
!
if(ltetra) then
wfshift = 0.5d0 * def(ipert) * dfpt_tetra_delta(ibnd,ik)
else
wfshift = 0.5d0 * def(ipert) * w0gauss( (ef-et(ibnd,ik))/degauss, ngauss) / degauss
end if
!
IF (noncolin) THEN
call zaxpy (npwx*npol,wfshift,evc(1,ibnd),1,dpsi(1,ibnd),1)
ELSE
call zaxpy (npw, wfshift, evc(1,ibnd), 1, dpsi(1,ibnd), 1)
ENDIF
enddo
!
! writes corrected delta_psi to iunit iudwf, k=kpoint,
!
IF (nksq > 1 .OR. npert > 1) CALL save_buffer(dpsi, lrdwf, iudwf, nrec)
enddo
enddo
!
END SUBROUTINE ef_shift
do ipert = 1, npert
do is = 1, nspin_mag
call zaxpy (dffts%nnr, def(ipert), ldoss(1,is), 1, drhoscf(1,is,ipert), 1)
enddo
enddo
!
CALL stop_clock ('ef_shift_wfc')
!
END SUBROUTINE ef_shift_wfc
!-------------------------------------------------------------------------
!
END MODULE efermi_shift
!-------------------------------------------------------------------------

1
PHonon/PH/print_clock_ph.f90
View File

@ -108,6 +108,7 @@ subroutine print_clock_ph
call print_clock ('dv_of_drho')
call print_clock ('mix_pot')
call print_clock ('ef_shift')
call print_clock ('ef_shift_wfc')
call print_clock ('localdos')
#if defined(__MPI)
call print_clock ('psymdvscf')

19
PHonon/PH/solve_linter.f90
View File

@ -73,7 +73,7 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
USE mp_pools, ONLY : inter_pool_comm
USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp
USE mp, ONLY : mp_sum
USE efermi_shift, ONLY : ef_shift, def
USE efermi_shift, ONLY : ef_shift, ef_shift_wfc, def
USE lrus, ONLY : int3_paw, becp1, int3_nc
USE lr_symm_base, ONLY : irotmq, minus_q, nsymq, rtau
USE eqv, ONLY : dvpsi
@ -401,14 +401,14 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
!
IF (okpaw) THEN
IF (lmetq0) CALL ef_shift(.FALSE., npe, dos_ef, ldos, ldoss, drhoscfh, &
IF (lmetq0) CALL ef_shift(npe, dos_ef, ldos, drhoscfh, &
dbecsum, becsum1, irr, sym_def)
DO ipert=1,npe
dbecsum(:,:,:,ipert)=2.0_DP *dbecsum(:,:,:,ipert) &
+becsumort(:,:,:,imode0+ipert)
ENDDO
ELSE
IF (lmetq0) CALL ef_shift(.FALSE., npe, dos_ef, ldos, ldoss, drhoscfh, &
IF (lmetq0) CALL ef_shift(npe, dos_ef, ldos, drhoscfh, &
irr=irr, sym_def=sym_def)
ENDIF
!
@ -453,9 +453,9 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
! And we mix with the old potential
!
IF (okpaw) THEN
!
! In this case we mix also dbecsum
!
!
! In this case we mix also dbecsum
!
call setmixout(npe*dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag*npe)/2, &
mixout, dvscfout, dbecsum, ndim, -1 )
call mix_potential (2*npe*dfftp%nnr*nspin_mag+2*ndim, mixout, mixin, &
@ -463,15 +463,14 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
nmix_ph, flmixdpot, convt)
call setmixout(npe*dfftp%nnr*nspin_mag,(nhm*(nhm+1)*nat*nspin_mag*npe)/2, &
mixin, dvscfin, dbecsum, ndim, 1 )
IF (lmetq0 .AND. convt) CALL ef_shift(.TRUE., npe, dos_ef, ldos, ldoss, drhoscf, &
dbecsum, becsum1, irr, sym_def)
ELSE
call mix_potential (2*npe*dfftp%nnr*nspin_mag, dvscfout, dvscfin, &
alpha_mix(kter), dr2, npe*tr2_ph/npol, iter, &
nmix_ph, flmixdpot, convt)
IF (lmetq0 .AND. convt) CALL ef_shift(.TRUE., npe, dos_ef, ldos, ldoss, drhoscf, &
irr=irr, sym_def=sym_def)
ENDIF
!
IF (lmetq0 .AND. convt) CALL ef_shift_wfc(npe, ldoss, drhoscf)
!
! check that convergent have been reached on ALL processors in this image
CALL check_all_convt(convt)