quantum-espresso/PHonon/PH/compute_weight.f90

124 lines
3.6 KiB
Fortran

!
! Copyright (C) 2001-2018 Quantum ESPRESSO 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 compute_weight (wgg)
!-----------------------------------------------------------------------
!! This routine implements Eq.(B19) of PRB 64, 235118 (2001). It computes the
!! weight to give to the v,v' terms in the orthogonality term.
!! The weights are computed for each k point.
!
USE kinds, ONLY : DP
USE klist, ONLY : wk, lgauss, degauss, ngauss, ltetra, degauss_cond
USE ener, ONLY : ef, ef_cond
USE wvfct, ONLY : nbnd, wg, et, nbnd_cond
USE paw_variables, ONLY : okpaw
USE qpoint, ONLY : nksq, ikks, ikqs
USE control_lr, ONLY : rec_code_read
USE dfpt_tetra_mod, ONLY : dfpt_tetra_ttheta
USE ldaU, ONLY : lda_plus_u
USE two_chem, ONLY : twochem
!
implicit none
!
real(DP) :: wgg(nbnd,nbnd,nksq)
!! output: the weights
!
! ... local variables
!
integer :: ik, ikk, ikq, ibnd, jbnd
! counters
real(DP) :: wg1, wg2, theta
! auxiliary variables
real(DP), external :: wgauss
real(DP), parameter :: eps = 1.0d-12
!
! DFPT+U: in the recover case the weights have to be computed in
! the routine dyn_hub_scf routine.
!
IF (lda_plus_u) GO TO 100
!
IF ((rec_code_read >= -20) .AND. (.NOT.okpaw)) RETURN
!
100 CONTINUE
!
IF(ltetra) THEN
!
DO ik = 1, nksq
DO ibnd = 1,nbnd
DO jbnd = 1, nbnd
wgg(ibnd, jbnd, ik) = dfpt_tetra_ttheta(jbnd, ibnd, ikks(ik))
END DO
END DO
END DO
!
return
!
END IF
!
do ik = 1, nksq
ikk = ikks(ik)
ikq = ikqs(ik)
!
! each band v ...
!
do ibnd = 1, nbnd
if (wk (ikk) .eq.0.d0) then
wg1 = 0.d0
else
wg1 = wg (ibnd, ikk) / wk (ikk)
!two chemical potentials case
if (twochem) then
if (ibnd.gt.(nbnd-nbnd_cond)) then
wg1 = wgauss ( (ef_cond - et (ibnd, ikk) ) / degauss_cond, ngauss)
else
wg1 = wgauss ( (ef - et (ibnd, ikk) ) / degauss, ngauss)
end if
end if
endif
!
! and each band v' ...
!
do jbnd = 1, nbnd
if (lgauss) then
theta = wgauss ( (et (jbnd,ikq) - et (ibnd,ikk) ) / degauss, 0)
wg2 = wgauss ( (ef - et (jbnd, ikq) ) / degauss, ngauss)
! two chemical potentials case
if (twochem) then
if (jbnd.gt.(nbnd-nbnd_cond)) then
wg2 = wgauss ( (ef_cond - et (jbnd, ikq) ) / degauss_cond, ngauss)
end if
end if
else
IF (et (jbnd,ikq) > et (ibnd,ikk)) THEN
theta = 1.0d0
ELSE
theta = 0.d0
ENDIF
IF (ABS(et (jbnd,ikq) - et (ibnd,ikk)) < 1.d-8) theta=0.5d0
if (wk (ikk) .le.eps) then
wg2 = 0.d0
else
wg2 = wg (jbnd, ikk) / wk (ikk)
endif
endif
wgg (ibnd, jbnd, ik) = wg1 * (1.d0 - theta) + wg2 * theta
enddo
enddo
! do ibnd=1,nbnd
! do jbnd=1,nbnd
! WRITE( stdout,'(3i5,f20.10)') ibnd, jbnd, ik,wgg(ibnd,jbnd,ik)
! enddo
! enddo
enddo
! call stop_ph(.true.)
return
end subroutine compute_weight