! Copyright (C) 2004-2009 Andrea Benassi and 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 .
!------------------------------
MODULE grid_module
!------------------------------
USE kinds, ONLY : DP
IMPLICIT NONE
PRIVATE
!
! general purpose vars
!
REAL(DP), ALLOCATABLE :: focc(:,:), wgrid(:)
REAL(DP) :: alpha
!
!
PUBLIC :: grid_build, grid_destroy
PUBLIC :: focc, wgrid, alpha
!
CONTAINS
!---------------------------------------------
SUBROUTINE grid_build(nw, wmax, wmin)
!-------------------------------------------
!
USE kinds, ONLY : DP
USE wvfct, ONLY : nbnd, wg
USE klist, ONLY : nks, wk, nelec, xk
USE lsda_mod, ONLY : nspin
USE uspp, ONLY : okvan
USE io_global, ONLY : stdout, ionode, ionode_id
!
IMPLICIT NONE
!
! input vars
INTEGER, INTENT(in) :: nw
REAL(DP), INTENT(in) :: wmax ,wmin
!
! local vars
INTEGER :: iw,ik,i,ierr
!
! check on the number of bands: we need to include empty bands in order to allow
! to write the transitions
!
!print*, nelec, nbnd
IF ( REAL(nbnd, DP) <= nelec / 2.0_DP ) CALL errore('epsilon', 'bad band number', 1)
!
! spin is not implemented
!
IF( nspin > 2 ) CALL errore('grid_build','Non collinear spin calculation not implemented',1)
!
! USPP are not implemented (dipole matrix elements are not trivial at all)
!
IF ( okvan ) CALL errore('grid_build','USPP are not implemented',1)
ALLOCATE ( focc( nbnd, nks), STAT=ierr )
IF (ierr/=0) CALL errore('grid_build','allocating focc', abs(ierr))
!
ALLOCATE( wgrid( nw ), STAT=ierr )
IF (ierr/=0) CALL errore('grid_build','allocating wgrid', abs(ierr))
!
! check on k point weights, no symmetry operations are allowed
!
! DO ik = 2, nks
! !
! IF ( abs( wk(1) - wk(ik) ) > 1.0d-8 ) &
! CALL errore('grid_build','non unifrom kpt grid', ik )
! !
! ENDDO
!
! PRINT the k point grid and weights as a test
!
!DO ik = 1, nks
! !
! IF ( ionode ) WRITE (6,*) xk(1,ik), xk(2,ik),xk(3,ik) , wk(ik)
! !
!ENDDO
!
! occupation numbers, to be normalized differently
! whether we are spin resolved or not
!
IF(nspin==1) THEN
DO ik = 1,nks
DO i = 1,nbnd
focc(i,ik)= wg(i, ik ) * 2.0_DP / wk( ik )
ENDDO
ENDDO
ELSEIF(nspin==2) THEN
DO ik = 1,nks
DO i = 1,nbnd
focc(i,ik)= wg(i, ik ) * 1.0_DP / wk( ik )
ENDDO
ENDDO
ENDIF
!
! set the energy grid
!
alpha = (wmax - wmin) / REAL(nw, DP)
!
DO iw = 1, nw
wgrid(iw) = wmin + iw * alpha
ENDDO
!
END SUBROUTINE grid_build
!
!
!----------------------------------
SUBROUTINE grid_destroy
!----------------------------------
IMPLICIT NONE
INTEGER :: ierr
!
IF ( allocated( focc) ) THEN
!
DEALLOCATE ( focc, wgrid, STAT=ierr)
CALL errore('grid_destroy','deallocating grid stuff',abs(ierr))
!
ENDIF
!
END SUBROUTINE grid_destroy
END MODULE grid_module
!------------------------------
PROGRAM epsilon
!------------------------------
!
! Compute the complex macroscopic dielectric function,
! at the RPA level, neglecting local field effects.
! Eps is computed both on the real or immaginary axis
!
! Authors: Andrea Benassi, Andrea Ferretti, Carlo Cavazzoni
!
! NOTE: Part of the basic implementation is taken from pw2gw.f90;
!
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout, ionode, ionode_id
USE mp, ONLY : mp_bcast
USE iotk_module
USE xml_io_base
USE io_files, ONLY : tmp_dir, prefix, outdir
USE constants, ONLY : RYTOEV
USE ener, ONLY : ef
USE klist, ONLY : lgauss
USE ktetra, ONLY : ltetra
USE wvfct, ONLY : nbnd
USE lsda_mod, ONLY : nspin
USE mp_global, ONLY : mp_startup
USE environment, ONLY : environment_start
!
IMPLICIT NONE
!
CHARACTER(LEN=256), EXTERNAL :: trimcheck
!
! input variables
!
INTEGER :: nw,nbndmin,nbndmax,nshell,ibndmin,ibndmax
REAL(DP) :: intersmear,intrasmear,wmax,wmin,shift,eta, qmod_par
CHARACTER(10) :: calculation,smeartype
LOGICAL :: metalcalc
!
NAMELIST / inputpp / prefix, outdir, calculation
NAMELIST / energy_grid / smeartype,intersmear,intrasmear,wmax,wmin,nbndmin,nbndmax,nw,shift,nshell,eta,ibndmin,ibndmax, qmod_par
!
! local variables
!
INTEGER :: ios
!---------------------------------------------
! program body
!---------------------------------------------
!
! initialise environment
!
#ifdef __MPI
CALL mp_startup ( )
#endif
CALL environment_start ( 'epsilon' )
!
! Set default values for variables in namelist
!
calculation = 'eps'
prefix = 'pwscf'
shift = 0.0d0
outdir = './'
intersmear = 0.136
wmin = 0.0d0
wmax = 30.0d0
eta = 0.3
nbndmin = 1
ibndmin = 1
nshell = 0
nbndmax = 0
ibndmax = 0
qmod_par = 0.5
nw = 600
smeartype = 'gauss'
intrasmear = 0.0d0
metalcalc = .false.
!
! this routine allows the user to redirect the input using -input
! instead of <
!
CALL input_from_file( )
!
! read input file
!
IF (ionode) WRITE( stdout, "( 2/, 5x, 'Reading input file...' ) " )
ios = 0
!
IF ( ionode ) READ (5, inputpp, IOSTAT=ios)
!
CALL mp_bcast ( ios, ionode_id )
!
IF (ios/=0) CALL errore('epsilon', 'reading namelist INPUTPP', abs(ios))
!
IF ( ionode ) THEN
!
READ (5, energy_grid, IOSTAT=ios)
!
tmp_dir = trimcheck(outdir)
!
ENDIF
!
CALL mp_bcast ( ios, ionode_id )
IF (ios/=0) CALL errore('epsilon', 'reading namelist ENERGY_GRID', abs(ios))
!
! ... Broadcast variables
!
IF (ionode) WRITE( stdout, "( 5x, 'Broadcasting variables...' ) " )
CALL mp_bcast( smeartype, ionode_id )
CALL mp_bcast( calculation, ionode_id )
CALL mp_bcast( prefix, ionode_id )
CALL mp_bcast( tmp_dir, ionode_id )
CALL mp_bcast( shift, ionode_id )
CALL mp_bcast( outdir, ionode_id )
CALL mp_bcast( intrasmear, ionode_id )
CALL mp_bcast( intersmear, ionode_id)
CALL mp_bcast( wmax, ionode_id )
CALL mp_bcast( wmin, ionode_id )
CALL mp_bcast( nw, ionode_id )
CALL mp_bcast( ibndmin, ionode_id )
CALL mp_bcast( ibndmax, ionode_id )
CALL mp_bcast( qmod_par, ionode_id )
CALL mp_bcast( nbndmin, ionode_id )
CALL mp_bcast( nbndmax, ionode_id )
CALL mp_bcast( nshell, ionode_id )
CALL mp_bcast( eta, ionode_id )
!
! read PW simulation parameters from prefix.save/data-file.xml
!
IF (ionode) WRITE( stdout, "( 5x, 'Reading PW restart file...' ) " )
CALL read_file
CALL openfil_pp
!
! few conversions
!
IF (ionode) WRITE(stdout,"(2/, 5x, 'Fermi energy [eV] is: ',f8.5)") ef *RYTOEV
IF (lgauss .or. ltetra) THEN
metalcalc=.true.
IF (ionode) WRITE( stdout, "( 5x, 'The system is a metal...' ) " )
ELSE
IF (ionode) WRITE( stdout, "( 5x, 'The system is a dielectric...' ) " )
ENDIF
IF (nbndmax == 0) nbndmax = nbnd
IF (ibndmax == 0) ibndmax = nbnd
!
! ... run the specific pp calculation
!
IF (ionode) WRITE(stdout,"(/, 5x, 'Performing ',a,' calculation...')") trim(calculation)
CALL start_clock( 'calculation' )
!
SELECT CASE ( trim(calculation) )
!
CASE ( 'gwppa' )
!
CALL gwppa_calc ( intersmear,intrasmear,nw,wmax,wmin,nbndmin,nbndmax,shift,metalcalc,nspin,nshell,eta,ibndmin,ibndmax,qmod_par)
!
CASE ( 'eps' )
!
CALL eps_calc ( intersmear,intrasmear,nw,wmax,wmin,nbndmin,nbndmax,shift,metalcalc,nspin )
!
CASE ( 'jdos' )
!
CALL jdos_calc ( smeartype,intersmear,nw,wmax,wmin,nbndmin,nbndmax,shift,nspin )
!
CASE ( 'offdiag' )
!
CALL offdiag_calc ( intersmear,intrasmear,nw,wmax,wmin,nbndmin,nbndmax,shift,metalcalc,nspin )
!
CASE ( 'occ' )
!
CALL occ_calc ()
!
CASE DEFAULT
!
CALL errore('epsilon','invalid CALCULATION = '//trim(calculation),1)
!
END SELECT
!
CALL stop_clock( 'calculation' )
!
! few info about timing
!
CALL stop_clock( 'epsilon' )
!
IF ( ionode ) WRITE( stdout , "(/)" )
!
CALL print_clock( 'epsilon' )
CALL print_clock( 'calculation' )
CALL print_clock( 'dipole_calc' )
!
IF ( ionode ) WRITE( stdout, * )
!
!
CALL stop_pp ()
END PROGRAM epsilon
!-----------------------------------------------------------------------------
SUBROUTINE gwppa_calc ( intersmear,intrasmear, nw, wmax, wmin, nbndmin, nbndmax, shift, &
metalcalc , nspin, nshell,eta, ibndmin, ibndmax, qmod_par)
!-----------------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE constants, ONLY : PI, RYTOEV, e2
USE cell_base, ONLY : tpiba2, omega, at, alat
USE wvfct, ONLY : nbnd, et, npw, igk, npwx, g2kin, ecutwfc
USE ener, ONLY : efermi => ef
USE klist, ONLY : nks, nkstot, degauss,xk,wk
USE io_global, ONLY : ionode, stdout
!
USE grid_module, ONLY : alpha, focc, wgrid, grid_build, grid_destroy
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm,me_bgrp
USE mp, ONLY : mp_sum
USE scf, ONLY : rho, rho_core, rhog_core, scf_type, v
USE wavefunctions_module, ONLY : evc !,psic
USE fft_base, ONLY : dfftp, dffts
USE gvect, ONLY : ngm, g
USE gvecs, ONLY : nls, nlsm
USE fft_interfaces, ONLY : fwfft, invfft
USE io_files, ONLY : nwordwfc, iunwfc
USE funct, ONLY : xc
USE grid_subroutines, ONLY : realspace_grids_info
!
IMPLICIT NONE
!
! input variables
!
INTEGER, INTENT(in) :: nw,nbndmin,nbndmax,nspin,nshell,ibndmin,ibndmax
REAL(DP), INTENT(in) :: wmax, wmin, intersmear,intrasmear, shift, eta, qmod_par
LOGICAL, INTENT(in) :: metalcalc
!
! local variables
!
INTEGER :: i,j,k, ik, iband1, iband2,is, jbnd, ibnd,ig, jg, kg, ir
INTEGER :: iw, iwp, ierr, nw1, nksigma, sysd
REAL(DP) :: etrans, const, w, renorm(3), wplasma, wplasma0, wtilde0, eps0, eps1, q, kTF, qG, qG1
COMPLEX(DP) :: wtilde
REAL(DP) :: wmin1, wmax1, wksigma, epssum, epsinv, epssum1
REAL(DP) :: vtxc, etxc, GN_freq, epsqdep
COMPLEX(DP) :: vxc,vc1,ovlp
COMPLEX(DP) :: ZDOTC
REAL(DP) :: inv_nr1, inv_nr2, inv_nr3
COMPLEX(DP) :: prefactor
REAL(DP) :: r(3), g2(3)
!
REAL(DP), ALLOCATABLE :: epsr(:,:), epsi(:,:), epsrc(:,:,:), epsic(:,:,:)
REAL(DP), ALLOCATABLE :: ieps(:,:), eels(:,:), iepsc(:,:,:), eelsc(:,:,:), ieps_qdep(:,:)
REAL(DP), ALLOCATABLE :: ieps_tmp (:,:)
REAL(DP), ALLOCATABLE :: dipole(:,:,:)
COMPLEX(DP), ALLOCATABLE :: sigmare (:)
! COMPLEX(DP), sigmax
REAL(DP), ALLOCATABLE :: spectrum(:)
COMPLEX(DP),ALLOCATABLE :: dipole_aux(:,:,:)
COMPLEX(DP) :: tempevc(dffts%nnr), tempevc1(dffts%nnr), tempevc2(dffts%nnr)
COMPLEX(DP),ALLOCATABLE :: psi(:),vpsi(:) !, tempevc(:)
COMPLEX(DP) :: evc1(npwx,nbnd)
REAL(DP) :: rhox, arhox, zeta, amag, vs, ex, ec, vx(2), vc(2), rhoneg(2), qmin
REAL(DP), PARAMETER :: vanishing_charge = 1.D-10
INTEGER :: npw1, ncount, indexqmin, indexgmin, npar, ipar, ipole
real(DP), ALLOCATABLE :: par (:), par_tmp (:)
INTEGER :: igk1 (npwx)
REAL(DP) :: absg0vec,ovlpsum
INTEGER, ALLOCATABLE :: gmap (:,:)
REAL(DP), ALLOCATABLE :: g0vec (:,:)
REAL(DP) :: g0 (3), gtmp(3,ngm)
REAL(DP) :: eta_large(2)
INTEGER :: ip, index, index0, ir_end, ng0vec
CHARACTER*2 label
CHARACTER*12 filesigma
CHARACTER*15 filespectrum
!
!--------------------------
! main routine body
!--------------------------
eta_large(1) = 0.d0
eta_large(2) = 0.d0
CALL grid_build(nw, wmax, wmin)
IF (nspin == 1) THEN
! Read PPA paramenter from file (named ppa.in)
! npar = 2
open (22,file="ppa.in",status="old",action="read")
read(22,*) npar
ALLOCATE(par(npar))
ALLOCATE(par_tmp (npar))
do ipar = 1, npar
read(22,*) par(ipar)
enddo
close(22)
! par(2) = sqrt(par(2))
! par(4) = sqrt(par(4))
! par(1) = -par(1)
! par(3) = -par(3)
!print *, par
! !----------------------------------------------------------------------------------
!print out the k-point grid
if(ionode)then
open(44,file='kgrid.dat')
write(44,'(A44)') '--------------------------------------------------------------------------------------'
write(44,'(A4,A10,A10,A10,A10)') 'ik', 'xk(1,ik)','xk(2,ik)','xk(3,ik)','wk(ik)'
write(44,'(A44)') '--------------------------------------------------------------------------------------'
do ik =1 ,nks
write(44,'(I4,F10.6,F10.6,F10.6,F10.6)') ik, xk(1,ik), xk(2,ik), xk(3,ik),wk(ik)
enddo
close(44)
endif
! !----------------------------------------------------------------------------------
!construct the frequency grid ------------------------------------
CALL grid_destroy()
wmin1 = -40.0d0 ! eV
wmax1 = 30.0d0 ! eV
nw1 = 500
CALL grid_build(nw1, wmax1, wmin1)
!-----------------------------------------------------------------
! print*, wgrid
call print_epsilon (1, par, wgrid, nw1 )
ALLOCATE( vpsi ( npwx ), STAT=ierr )
IF (ierr/=0) CALL errore('gwppa','allocating vpsi', abs(ierr))
ALLOCATE( psi ( npwx ), STAT=ierr )
IF (ierr/=0) CALL errore('gwppa','allocating psi', abs(ierr))
ALLOCATE( sigmare (nw1), STAT=ierr )
IF (ierr/=0) CALL errore('gwppa','allocating sigmare', abs(ierr))
ALLOCATE( spectrum (nw1), STAT=ierr )
IF (ierr/=0) CALL errore('gwppa','allocating spectrum', abs(ierr))
nksigma = 1
wksigma = 2.0d0 / (nks - nksigma)
!-----------------------------------------------------------------
! EXTERNAL LOOP OVER BANDS FOR WHICH THE SELF-ENERGY MUST BE CALCULATED
! DO ibnd = 4,4
! IF(ionode) PRINT*, "Smallest band index: ", ibndmin
! IF(ionode) PRINT*, "tpiba2: ", tpiba2,alat
! IF(ionode) PRINT*, "Biggest band index: ", ibndmax
! IF(ionode) PRINT*, "Parameter for q-dependence: ", qmod_par
ALLOCATE( ieps_tmp(3,nw) )
DO ibnd = ibndmin, ibndmax
IF(ionode) PRINT*, "Evaluating sigma for band: ", ibnd
!calculate expectation value of v_xc ----------------------------
ik=1
CALL gk_sort (xk (1, ik), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
IF (ibnd .LT.10) THEN
WRITE(label,'(A,I1)') "0", ibnd
ELSE
WRITE(label,'(I2)') ibnd
ENDIF
filespectrum = "spectrum_"//label//".dat"
filesigma = "sigma_"//label//".dat"
tempevc (:) = (0.0d0,0.0d0)
vpsi (:) = (0.0d0,0.0d0)
vxc = (0.0d0,0.0d0)
vc1 = (0.0d0,0.0d0)
! !construct the exchange-correlation potential ------------------
v%of_r(:,:) = (0.0d0, 0.0d0)
tempevc(:) = (0.0d0, 0.0d0)
vpsi (:) = (0.0d0, 0.0d0)
CALL davcio (evc, nwordwfc, iunwfc, ik, -1)
DO ig = 1 , npw
tempevc(nls(igk(ig))) = evc(ig,ibnd)
ENDDO
CALL invfft ('Wave', tempevc(:), dffts)
!evaluate matrix elements of v_xc
!v%of_r(:,:) = (0.0d0, 0.0d0)
CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, v%of_r )
DO ir = 1, dfftp%nnr
tempevc (ir) = v%of_r (ir, 1) * tempevc (ir)
ENDDO
CALL fwfft ('Wave', tempevc(:), dffts)
DO ig = 1, npw
vpsi(ig) = tempevc (nls(igk(ig)))
ENDDO
vxc = ZDOTC (npw, evc (1, ibnd), 1, vpsi, 1)
CALL mp_sum( vxc, intra_pool_comm )
! if (ionode) PRINT * , 'VXC = ', real(vxc)*RYTOEV, 'VC = ', real(vc1)*RYTOEV
! !--------------------------------------------
!construct the exchange-correlation potential ------------------
v%of_r(:,:) = (0.0d0, 0.0d0)
tempevc(:) = (0.0d0, 0.0d0)
vpsi (:) = (0.0d0, 0.0d0)
CALL davcio (evc, nwordwfc, iunwfc, ik, -1)
DO ig = 1 , npw
tempevc(nls(igk(ig))) = evc(ig,ibnd)
ENDDO
CALL invfft ('Wave', tempevc(:), dffts)
!evaluate matrix elements of v_xc
!v%of_r(:,:) = (0.0d0, 0.0d0)
CALL get_vc( rho, rho_core, rhog_core, etxc, vtxc, v%of_r )
DO ir = 1, dfftp%nnr
tempevc (ir) = v%of_r (ir, 1) * tempevc (ir)
ENDDO
CALL fwfft ('Wave', tempevc(:), dffts)
DO ig = 1, npw
vpsi(ig) = tempevc (nls(igk(ig)))
ENDDO
vc1 = ZDOTC (npw, evc (1, ibnd), 1, vpsi, 1)
!CALL mp_sum( vxc, intra_pool_comm )
CALL mp_sum( vc1, intra_pool_comm )
if (ionode) PRINT * , 'VXC = ', real(vxc)*RYTOEV, 'VC = ', real(vc1)*RYTOEV
!--------------------------------------------
IF ( ionode ) WRITE(6,*) 'EVHERE:', et(ibnd,1) * RYTOEV, efermi*RYTOEV
!evaluate the super-approximated gwppa self-energy ----------------------------
sigmare (:) = (0.0d0,0.0d0)
vxc = (0.0d0,0.0d0)
!vc1 = (0.0d0,0.0d0)
!nshell = 1
ng0vec = (2*nshell+1)**3
IF(.NOT.ALLOCATED(gmap))THEN
ALLOCATE (gmap(ngm,ng0vec))
ENDIF
IF(.NOT.ALLOCATED(g0vec))THEN
ALLOCATE (g0vec(3,ng0vec))
ENDIF
gmap(:,:)=0.d0
CALL refold (gmap,nshell,g0vec)
!determine the minimum q-point, for the integration of the q=0 singularity
qmin = 100.d0
DO ik = nksigma+1, nks
DO ig = 1 , ng0vec
q =sum((xk(:,ik)-xk(:,1)+g0vec(:,ig))**2.d0)
IF(q .LT. qmin)THEN
qmin = q
indexqmin = ik
indexgmin = ig
ENDIF
ENDDO
ENDDO
!-----------------------------------------------------------------
CALL gk_sort (xk (1, 1), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio (evc, nwordwfc, iunwfc, 1, -1)
tempevc(:)=(0.d0,0.d0)
tempevc(nls(igk(1:npw))) = evc(1:npw,ibnd)
! !test gmap
! DO ig=1,1
! DO jg =1,(2*nshell+1)**3
! IF(ionode.and..not. gmap(ig,jg).eq.0)THEN
! PRINT*, '--------------------------'
! PRINT*, ig, jg ,igk(ig), gmap(ig,jg)
! PRINT*, g(:,ig)
! PRINT*, g(:,gmap(ig,jg))
! ENDIF
! ENDDO
! ENDDO
!DO ik = 1, nks !always exclude the gamma (q=0) point from the summation
DO ik = nksigma+1, nks !always exclude the gamma (q=0) point from the summation
! IF (ionode) PRINT*, ' ik = ' , ik , ' / ', nks
!construct the overlap --------------------------------------
CALL gk_sort (xk (1,ik), ngm, g, ecutwfc / tpiba2, npw1, igk1, g2kin)
CALL davcio (evc1, nwordwfc, iunwfc, ik, -1)
DO ig = 1 , ng0vec !loop over G-vectors
DO jbnd = 1, nbnd
!!! !!! !calculate an overlap by hand --------------------
!!! !!! ! IF(ik.EQ.1)THEN
!!! !!! ovlp = (0.d0,0.d0)
!!! !!! g0(:) = g0vec(:,ig)
!!! !!! gtmp(:,:) = g(:,:)
!!! !!!
!!! !!! call cryst_to_cart (ngm, gtmp, at, -1)
!!! !!! call cryst_to_cart (1, g0, at, -1)
!!! !!!
!!! !!! DO kg = 1, npw
!!! !!! g2(:) = gtmp(:,igk(kg))+g0(:)
!!! !!! DO jg = 1, npw1
!!! !!! IF( (nint(gtmp(1,igk1(jg))) .EQ. nint(g2(1))) .AND. &
!!! !!! (nint(gtmp(2,igk1(jg))) .EQ. nint(g2(2))) .AND. &
!!! !!! (nint(gtmp(3,igk1(jg))) .EQ. nint(g2(3))) )THEN
!!! !!! ovlp = ovlp + evc(kg,ibnd) * conjg(evc1(jg,jbnd))
!!! !!! ENDIF
!!! !!! ENDDO
!!! !!! ENDDO
!!! !!!
!!! !!! IF(ionode)PRINT*, 'ibnd = ', ibnd, 'jbnd = ', jbnd
!!! !!! IF(ionode)PRINT*, 'ovlp1 = ', ovlp!, ncount
!!! !!! ! ENDIF
!!! !!! !-------------------------------------------------
ovlp = (0.d0,0.d0)
tempevc1(:)=(0.d0,0.d0)
DO jg = 1, npw1 !loop over G'-vectors
IF(gmap(igk1(jg),ig) .GT. 0 )THEN
tempevc1(nls(gmap(igk1(jg),ig))) = evc1(jg,jbnd)
ENDIF
ENDDO
ovlp = ZDOTC (dffts%nnr, tempevc1, 1, tempevc, 1)
CALL mp_sum( ovlp, intra_pool_comm )
!------------------------------------------------------------
qG = sum((xk(:,ik)-xk(:,1)+g0vec(:,ig))**2) *tpiba2
qG1 = sum((xk(:,ik) +g0vec(:,ig))**2)
DO ipole = 1, npar/2
wtilde = (par( ipole * 2 ) + (0.d0,1.d0) * eta_large(ipole) ) * dsqrt(1.d0 + qG / 1.2 )
!if( ionode .and. ibnd == 1 .and. jbnd == 1 ) print *, qG, qG / tpiba2, wtilde
!wtilde = par( 2 + ( ipole - 1 ) * 2 ) * dsqrt(1.d0 + qG1 / qmod_par)
! if( ionode .and. ibnd == 1 .and. jbnd == 1 ) print * , tpiba2
! if( ionode .and. ibnd == 1 .and. jbnd == 1 ) print * , qG1, dsqrt(1.d0 + qG1 /0.6), ik
wplasma = par( ipole * 2 - 1 )
prefactor = ovlp * conjg (ovlp) * &
wplasma**2 / (2.0d0*wtilde ) * 4.0d0 * PI
IF (.NOT. (ik.EQ.indexqmin .AND. ig .EQ. indexgmin)) THEN
DO iw = 1, nw1
sigmare(iw) = sigmare(iw) + &
prefactor / (wgrid(iw) + SIGN(real(wtilde),efermi - et(jbnd,ik)) - et(jbnd,ik) &
* RYTOEV - (0.0d0,1.0d0) * SIGN(eta,efermi - et(jbnd,ik))) &
/ qG * wksigma / omega *RYTOEV !* 1.d0/ (1.d0 + qG /kTF/tpiba2)
ENDDO
ELSE
!IF (ionode) WRITE(6,*) "I am neglecting the q=0 singularity!"
DO iw = 1, nw1
sigmare(iw) = sigmare(iw) + &
prefactor / ( wgrid(iw) + SIGN(real(wtilde),efermi - et(jbnd,ik)) - et(jbnd,ik) &
* RYTOEV - (0.0d0,1.0d0) * SIGN(eta,efermi - et(jbnd,ik))) * &
( 3.d0 * wksigma / 4.d0 / PI / omega )**(1.d0/3.d0)/PI * RYTOEV
ENDDO !loop over frequency points
ENDIF ! (.NOT.(ik.EQ.indexqmin .AND. absg0vec .LT. 1e-3))
ENDDO
ENDDO !loop over all bands
ENDDO !loop over G-vectors shells
ENDDO !loop over k-points
CALL mp_sum( sigmare, inter_pool_comm )
!------------------------------------------------------------------------------
!evaluate the spectral function -----------------------------------------------
DO iw = 1, nw1
spectrum(iw) = abs(aimag (sigmare(iw))) / &
((wgrid(iw) - (et(ibnd,1)-real(vc1)) * RYTOEV &
- real(sigmare(iw)))**2 + aimag (sigmare(iw))**2)
ENDDO
!------------------------------------------------------------------------------
IF( ionode ) PRINT*, ' v_c [eV]: ', real(vc1) * RYTOEV
!print to file ----------------------------------------------------------------
IF ( ionode ) THEN
OPEN (30, FILE=filesigma)
OPEN (31, FILE=filespectrum)
DO iw = 1, nw1
WRITE(30,*) wgrid(iw), real(sigmare(iw)), aimag (sigmare(iw))
WRITE(31,*) wgrid(iw)-efermi*RYTOEV, spectrum(iw)
!WRITE(31,*) wgrid(iw)-efermi*RYTOEV, spectrum(iw)
ENDDO
CLOSE(30)
CLOSE(31)
ENDIF
! DEALLOCATE (sigmare)
!------------------------------------------------------------------------------
ENDDO ! LOOP OVER OCCUPIED BANDS
ELSEIF (nspin == 2 ) THEN
IF ( ionode ) WRITE(6,*) 'GWPPA works only for nspin = 1 !'
ENDIF
IF(ionode ) WRITE(6,*) " DONE !!! "
!--------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------
!
! local cleaning
!
CALL grid_destroy()
!
DEALLOCATE ( dipole, dipole_aux )
END SUBROUTINE gwppa_calc
!-----------------------------------------------------------------------------
SUBROUTINE eps_calc ( intersmear,intrasmear, nw, wmax, wmin, nbndmin, nbndmax, shift, &
metalcalc , nspin)
!-----------------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE constants, ONLY : PI, RYTOEV
USE cell_base, ONLY : tpiba2, omega
USE wvfct, ONLY : nbnd, et
USE ener, ONLY : efermi => ef
USE klist, ONLY : nks, nkstot, degauss
USE io_global, ONLY : ionode, stdout
!
USE grid_module, ONLY : alpha, focc, wgrid, grid_build, grid_destroy
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm
USE mp, ONLY : mp_sum
!
IMPLICIT NONE
!
! input variables
!
INTEGER, INTENT(in) :: nw,nbndmin,nbndmax,nspin
REAL(DP), INTENT(in) :: wmax, wmin, intersmear,intrasmear, shift
LOGICAL, INTENT(in) :: metalcalc
!
! local variables
!
INTEGER :: i, ik, iband1, iband2,is
INTEGER :: iw, iwp, ierr
REAL(DP) :: etrans, const, w, renorm(3)
!
REAL(DP), ALLOCATABLE :: epsr(:,:), epsi(:,:), epsrc(:,:,:), epsic(:,:,:)
REAL(DP), ALLOCATABLE :: ieps(:,:), eels(:,:), iepsc(:,:,:), eelsc(:,:,:)
REAL(DP), ALLOCATABLE :: dipole(:,:,:)
COMPLEX(DP),ALLOCATABLE :: dipole_aux(:,:,:)
!
!--------------------------
! main routine body
!--------------------------
!
!
! perform some consistency checks, calculate occupation numbers and setup w grid
!
CALL grid_build(nw, wmax, wmin)
!
! allocate main spectral and auxiliary quantities
!
ALLOCATE( dipole(3, nbnd, nbnd), STAT=ierr )
IF (ierr/=0) CALL errore('epsilon','allocating dipole', abs(ierr) )
!
ALLOCATE( dipole_aux(3, nbnd, nbnd), STAT=ierr )
IF (ierr/=0) CALL errore('epsilon','allocating dipole_aux', abs(ierr) )
!
! spin unresolved calculation
!
IF (nspin == 1) THEN
!
ALLOCATE( epsr( 3, nw), epsi( 3, nw), eels( 3, nw), ieps(3,nw ), STAT=ierr )
IF (ierr/=0) CALL errore('epsilon','allocating eps', abs(ierr))
!
! initialize response functions
!
epsr(:,:) = 0.0_DP
epsi(:,:) = 0.0_DP
ieps(:,:) = 0.0_DP
!
! main kpt loop
!
kpt_loop: &
DO ik = 1, nks
!
! For every single k-point: order k+G for
! read and distribute wavefunctions
! compute dipole matrix 3 x nbnd x nbnd parallel over g
! recover g parallelism getting the total dipole matrix
!
CALL dipole_calc( ik, dipole_aux, metalcalc , nbndmin, nbndmax)
!
dipole(:,:,:)= tpiba2 * REAL( dipole_aux(:,:,:) * conjg(dipole_aux(:,:,:)), DP )
!
! Calculation of real and immaginary parts
! of the macroscopic dielettric function from dipole
! approximation.
! 'intersmear' is the brodening parameter
!
!Interband
!
DO iband2 = nbndmin,nbndmax
!
IF ( focc(iband2,ik) < 2.0d0) THEN
DO iband1 = nbndmin,nbndmax
!
IF (iband1==iband2) CYCLE
IF ( focc(iband1,ik) >= 1e-4 ) THEN
IF (abs(focc(iband2,ik)-focc(iband1,ik))< 1e-3) CYCLE
!
! transition energy
!
etrans = ( et(iband2,ik) -et(iband1,ik) ) * RYTOEV + shift
!
! loop over frequencies
!
DO iw = 1, nw
!
w = wgrid(iw)
!
epsi(:,iw) = epsi(:,iw) + dipole(:,iband1,iband2) * intersmear * w* &
RYTOEV**3 * (focc(iband1,ik))/ &
(( (etrans**2 -w**2 )**2 + intersmear**2 * w**2 )* etrans )
epsr(:,iw) = epsr(:,iw) + dipole(:,iband1,iband2) * RYTOEV**3 * &
(focc(iband1,ik)) * &
(etrans**2 - w**2 ) / &
(( (etrans**2 -w**2 )**2 + intersmear**2 * w**2 )* etrans )
ENDDO
ENDIF
ENDDO
ENDIF
ENDDO
!
!Intraband (only if metalcalc is true)
!
IF (metalcalc) THEN
DO iband1 = nbndmin,nbndmax
!
IF ( focc(iband1,ik) < 2.0d0) THEN
IF ( focc(iband1,ik) >= 1e-4 ) THEN
!
! loop over frequencies
!
DO iw = 1, nw
!
w = wgrid(iw)
!
epsi(:,iw) = epsi(:,iw) + dipole(:,iband1,iband1) * intrasmear * w* &
RYTOEV**2 * (exp((et(iband1,ik)-efermi)/degauss ))/ &
(( w**4 + intrasmear**2 * w**2 )*(1+exp((et(iband1,ik)-efermi)/ &
degauss))**2*degauss )
epsr(:,iw) = epsr(:,iw) - dipole(:,iband1,iband1) * RYTOEV**2 * &
(exp((et(iband1,ik)-efermi)/degauss )) * w**2 / &
(( w**4 + intrasmear**2 * w**2 )*(1+exp((et(iband1,ik)-efermi)/ &
degauss))**2*degauss )
ENDDO
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO kpt_loop
!
! recover over kpt parallelization (inter_pool)
!
CALL mp_sum( epsr, inter_pool_comm )
CALL mp_sum( epsi, inter_pool_comm )
!
! impose the correct normalization
!
const = 64.0d0 * PI / ( omega * REAL(nkstot, DP) )
epsr(:,:) = 1.0_DP + epsr(:,:) * const
epsi(:,:) = epsi(:,:) * const
!
! Calculation of eels spectrum
!
DO iw = 1, nw
!
eels(:,iw) = epsi(:,iw) / ( epsr(:,iw)**2 + epsi(:,iw)**2 )
!
ENDDO
!
! calculation of dielectric function on the immaginary frequency axe
!
DO iw = 1, nw
DO iwp = 2, nw
!
ieps(:,iw) = ieps(:,iw) + wgrid(iwp) * epsi(:,iwp) / ( wgrid(iwp)**2 + wgrid(iw)**2)
!
ENDDO
ENDDO
ieps(:,:) = 1.0d0 + 2 / PI * ieps(:,:) * alpha
!
! check dielectric function normalizzation via sumrule
!
DO i=1,3
renorm(i) = alpha * sum( epsi(i,:) * wgrid(:) )
ENDDO
!
IF ( ionode ) THEN
!
WRITE(stdout,"(/,5x, 'The bulk xx plasmon frequency [eV] is: ',f15.9 )") sqrt(renorm(1) * 2.0d0 / PI)
WRITE(stdout,"(5x, 'The bulk yy plasmon frequency [eV] is: ',f15.9 )") sqrt(renorm(2) * 2.0d0 / PI)
WRITE(stdout,"(5x, 'The bulk zz plasmon frequency [eV] is: ',f15.9 )") sqrt(renorm(3) * 2.0d0 / PI)
WRITE(stdout,"(/,5x, 'Writing output on file...' )")
!
! write results on data files
!
OPEN (30, FILE='epsr.dat', FORM='FORMATTED' )
OPEN (40, FILE='epsi.dat', FORM='FORMATTED' )
OPEN (41, FILE='eels.dat', FORM='FORMATTED' )
OPEN (42, FILE='ieps.dat', FORM='FORMATTED' )
!
WRITE(30, "(2x,'# energy grid [eV] epsr_x epsr_y epsr_z')" )
WRITE(40, "(2x,'# energy grid [eV] epsi_x epsi_y epsi_z')" )
WRITE(41, "(2x,'# energy grid [eV] eels components [arbitrary units]')" )
WRITE(42, "(2x,'# energy grid [eV] ieps_x ieps_y ieps_z ')" )
!
DO iw =1, nw
!
WRITE(30,"(4f15.6)") wgrid(iw), epsr(1:3, iw)
WRITE(40,"(4f15.6)") wgrid(iw), epsi(1:3, iw)
WRITE(41,"(4f15.6)") wgrid(iw), eels(1:3, iw)
WRITE(42,"(4f15.6)") wgrid(iw), ieps(1:3, iw)
!
ENDDO
!
CLOSE(30)
CLOSE(40)
CLOSE(41)
CLOSE(42)
!
ENDIF
DEALLOCATE ( epsr, epsi, eels, ieps)
!
! collinear spin calculation
!
ELSEIF (nspin == 2 ) THEN
!
ALLOCATE( epsrc( 0:1, 3, nw), epsic( 0:1,3, nw), eelsc( 0:1,3, nw), iepsc(0:1,3,nw ), STAT=ierr )
IF (ierr/=0) CALL errore('epsilon','allocating eps', abs(ierr))
!
! initialize response functions
!
epsrc(:,:,:) = 0.0_DP
epsic(:,:,:) = 0.0_DP
iepsc(:,:,:) = 0.0_DP
!
! main kpt loop
!
spin_loop: &
DO is=0,1
kpt_loopspin: &
! if nspin=2 the number of nks must be even (even if the calculation
! is performed at gamma point only), so nks must be always a multiple of 2
DO ik = 1 + is * int(nks/2), int(nks/2) + is * int(nks/2)
!
! For every single k-point: order k+G for
! read and distribute wavefunctions
! compute dipole matrix 3 x nbnd x nbnd parallel over g
! recover g parallelism getting the total dipole matrix
!
CALL dipole_calc( ik, dipole_aux, metalcalc , nbndmin, nbndmax)
!
dipole(:,:,:)= tpiba2 * REAL( dipole_aux(:,:,:) * conjg(dipole_aux(:,:,:)), DP )
!
! Calculation of real and immaginary parts
! of the macroscopic dielettric function from dipole
! approximation.
! 'intersmear' is the brodening parameter
!
!Interband
!
DO iband2 = nbndmin,nbndmax
!
IF ( focc(iband2,ik) < 1.0d0) THEN
DO iband1 = nbndmin,nbndmax
!
IF (iband1==iband2) CYCLE
IF ( focc(iband1,ik) >= 1e-4 ) THEN
IF (abs(focc(iband2,ik)-focc(iband1,ik))< 1e-3) CYCLE
!
! transition energy
!
etrans = ( et(iband2,ik) -et(iband1,ik) ) * RYTOEV + shift
!
! loop over frequencies
!
DO iw = 1, nw
!
w = wgrid(iw)
!
epsic(is,:,iw) = epsic(is,:,iw) + dipole(:,iband1,iband2) * intersmear * w* &
RYTOEV**3 * (focc(iband1,ik))/ &
(( (etrans**2 -w**2 )**2 + intersmear**2 * w**2 )* etrans )
epsrc(is,:,iw) = epsrc(is,:,iw) + dipole(:,iband1,iband2) * RYTOEV**3 * &
(focc(iband1,ik)) * &
(etrans**2 - w**2 ) / &
(( (etrans**2 -w**2 )**2 + intersmear**2 * w**2 )* etrans )
ENDDO
ENDIF
ENDDO
ENDIF
ENDDO
!
!Intraband (only if metalcalc is true)
!
IF (metalcalc) THEN
DO iband1 = nbndmin,nbndmax
!
IF ( focc(iband1,ik) < 1.0d0) THEN
IF ( focc(iband1,ik) >= 1e-4 ) THEN
!
! loop over frequencies
!
DO iw = 1, nw
!
w = wgrid(iw)
!
epsic(is,:,iw) = epsic(is,:,iw) + dipole(:,iband1,iband1) * intrasmear * w* &
RYTOEV**2 * (exp((et(iband1,ik)-efermi)/degauss ))/ &
(( w**4 + intrasmear**2 * w**2 )*(1+exp((et(iband1,ik)-efermi)/ &
degauss))**2*degauss )
epsrc(is,:,iw) = epsrc(is,:,iw) - dipole(:,iband1,iband1) * RYTOEV**2 * &
(exp((et(iband1,ik)-efermi)/degauss )) * w**2 / &
(( w**4 + intrasmear**2 * w**2 )*(1+exp((et(iband1,ik)-efermi)/ &
degauss))**2*degauss )
ENDDO
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO kpt_loopspin
ENDDO spin_loop
!
! recover over kpt parallelization (inter_pool)
!
CALL mp_sum( epsr, inter_pool_comm )
CALL mp_sum( epsi, inter_pool_comm )
!
! impose the correct normalization
!
const = 128.0d0 * PI / ( omega * REAL(nkstot, DP) )
epsrc(:,:,:) = 1.0_DP + epsrc(:,:,:) * const
epsic(:,:,:) = epsic(:,:,:) * const
!
! Calculation of eels spectrum
!
DO iw = 1, nw
!
eelsc(:,:,iw) = epsic(:,:,iw) / ( epsrc(:,:,iw)**2 + epsic(:,:,iw)**2 )
!
ENDDO
!
! calculation of dielectric function on the immaginary frequency axe
!
DO iw = 1, nw
DO iwp = 2, nw
!
iepsc(:,:,iw) = iepsc(:,:,iw) + wgrid(iwp) * epsic(:,:,iwp) / ( wgrid(iwp)**2 + wgrid(iw)**2)
!
ENDDO
ENDDO
iepsc(:,:,:) = 1.0d0 + 2.0_DP / PI * iepsc(:,:,:) * alpha
IF (ionode) THEN
WRITE(stdout,"(/,5x, 'Writing output on file...' )")
!
! write results on data files
!
OPEN (30, FILE='uepsr.dat', FORM='FORMATTED' )
OPEN (40, FILE='uepsi.dat', FORM='FORMATTED' )
OPEN (41, FILE='ueels.dat', FORM='FORMATTED' )
OPEN (42, FILE='uieps.dat', FORM='FORMATTED' )
OPEN (43, FILE='depsr.dat', FORM='FORMATTED' )
OPEN (44, FILE='depsi.dat', FORM='FORMATTED' )
OPEN (45, FILE='deels.dat', FORM='FORMATTED' )
OPEN (46, FILE='dieps.dat', FORM='FORMATTED' )
OPEN (47, FILE='epsr.dat', FORM='FORMATTED' )
OPEN (48, FILE='epsi.dat', FORM='FORMATTED' )
OPEN (49, FILE='eels.dat', FORM='FORMATTED' )
OPEN (50, FILE='ieps.dat', FORM='FORMATTED' )
!
WRITE(30, "(2x,'# energy grid [eV] epsr_x epsr_y epsr_z')" )
WRITE(40, "(2x,'# energy grid [eV] epsi_x epsi_y epsi_z')" )
WRITE(41, "(2x,'# energy grid [eV] eels components [arbitrary units]')" )
WRITE(42, "(2x,'# energy grid [eV] ieps_x ieps_y ieps_z ')" )
WRITE(43, "(2x,'# energy grid [eV] epsr_x epsr_y epsr_z')" )
WRITE(44, "(2x,'# energy grid [eV] epsi_x epsi_y epsi_z')" )
WRITE(45, "(2x,'# energy grid [eV] eels components [arbitrary units]')" )
WRITE(46, "(2x,'# energy grid [eV] ieps_x ieps_y ieps_z ')" )
WRITE(47, "(2x,'# energy grid [eV] epsr_x epsr_y epsr_z')" )
WRITE(48, "(2x,'# energy grid [eV] epsi_x epsi_y epsi_z')" )
WRITE(49, "(2x,'# energy grid [eV] eels components [arbitrary units]')" )
WRITE(50, "(2x,'# energy grid [eV] ieps_x ieps_y ieps_z ')" )
!
DO iw =1, nw
!
WRITE(30,"(4f15.6)") wgrid(iw), epsrc(0,1:3, iw)
WRITE(40,"(4f15.6)") wgrid(iw), epsic(0,1:3, iw)
WRITE(41,"(4f15.6)") wgrid(iw), eelsc(0,1:3, iw)
WRITE(42,"(4f15.6)") wgrid(iw), iepsc(0,1:3, iw)
WRITE(43,"(4f15.6)") wgrid(iw), epsrc(1,1:3, iw)
WRITE(44,"(4f15.6)") wgrid(iw), epsic(1,1:3, iw)
WRITE(45,"(4f15.6)") wgrid(iw), eelsc(1,1:3, iw)
WRITE(46,"(4f15.6)") wgrid(iw), iepsc(1,1:3, iw)
WRITE(47,"(4f15.6)") wgrid(iw), epsrc(1,1:3, iw)+epsrc(0,1:3, iw)
WRITE(48,"(4f15.6)") wgrid(iw), epsic(1,1:3, iw)+epsic(0,1:3, iw)
WRITE(49,"(4f15.6)") wgrid(iw), eelsc(1,1:3, iw)+eelsc(0,1:3, iw)
WRITE(50,"(4f15.6)") wgrid(iw), iepsc(1,1:3, iw)+iepsc(0,1:3, iw)
!
ENDDO
!
CLOSE(30)
CLOSE(40)
CLOSE(41)
CLOSE(42)
CLOSE(43)
CLOSE(44)
CLOSE(45)
CLOSE(46)
CLOSE(47)
CLOSE(48)
CLOSE(49)
CLOSE(50)
!
ENDIF
DEALLOCATE ( epsrc, epsic, eelsc, iepsc)
ENDIF
!
! local cleaning
!
CALL grid_destroy()
!
DEALLOCATE ( dipole, dipole_aux )
END SUBROUTINE eps_calc
!----------------------------------------------------------------------------------------
SUBROUTINE jdos_calc ( smeartype,intersmear,nw,wmax,wmin,nbndmin,nbndmax,shift,nspin )
!--------------------------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE constants, ONLY : PI, RYTOEV
USE wvfct, ONLY : nbnd, et
USE klist, ONLY : nks
USE io_global, ONLY : ionode, stdout
USE grid_module, ONLY : alpha, focc, wgrid, grid_build, grid_destroy
!
IMPLICIT NONE
!
! input variables
!
INTEGER, INTENT(in) :: nw,nbndmin,nbndmax,nspin
REAL(DP), INTENT(in) :: wmax, wmin, intersmear, shift
CHARACTER(*), INTENT(in) :: smeartype
!
! local variables
!
INTEGER :: ik, is, iband1, iband2
INTEGER :: iw, ierr
REAL(DP) :: etrans, w, renorm, count, srcount(0:1), renormzero,renormuno
!
REAL(DP), ALLOCATABLE :: jdos(:),srjdos(:,:)
!
!--------------------------
! main routine body
!--------------------------
!
! No wavefunctions are needed in order to compute jdos, only eigenvalues,
! they are distributed to each task so
! no mpi calls are necessary in this routine
!
! perform some consistency checks, calculate occupation numbers and setup w grid
!
CALL grid_build(nw, wmax, wmin )
!
! spin unresolved calculation
!
IF (nspin == 1) THEN
!
! allocate main spectral and auxiliary quantities
!
ALLOCATE( jdos(nw), STAT=ierr )
IF (ierr/=0) CALL errore('epsilon','allocating jdos',abs(ierr))
!
! initialize jdos
!
jdos(:)=0.0_DP
! Initialising a counter for the number of transition
count=0.0_DP
!
! main kpt loop
!
IF (smeartype=='lorentz') THEN
kpt_lor: &
DO ik = 1, nks
!
! Calculation of joint density of states
! 'intersmear' is the brodening parameter
!
DO iband2 = 1,nbnd
IF ( focc(iband2,ik) < 2.0d0) THEN
DO iband1 = 1,nbnd
!
IF ( focc(iband1,ik) >= 1.0d-4 ) THEN
!
! transition energy
!
etrans = ( et(iband2,ik) -et(iband1,ik) ) * RYTOEV + shift
!
IF( etrans < 1.0d-10 ) CYCLE
count = count + (focc(iband1,ik)-focc(iband2,ik))
!
! loop over frequencies
!
DO iw = 1, nw
!
w = wgrid(iw)
!
jdos(iw) = jdos(iw) + intersmear * (focc(iband1,ik)-focc(iband2,ik)) &
/ ( PI * ( (etrans -w )**2 + (intersmear)**2 ) )
ENDDO
ENDIF
ENDDO
ENDIF
ENDDO
ENDDO kpt_lor
ELSEIF (smeartype=='gauss') THEN
kpt_gauss: &
DO ik = 1, nks
!
! Calculation of joint density of states
! 'intersmear' is the brodening parameter
!
DO iband2 = 1,nbnd
DO iband1 = 1,nbnd
!
IF ( focc(iband2,ik) < 2.0d0) THEN
IF ( focc(iband1,ik) >= 1.0d-4 ) THEN
!
! transition energy
!
etrans = ( et(iband2,ik) -et(iband1,ik) ) * RYTOEV + shift
!
IF( etrans < 1.0d-10 ) CYCLE
! loop over frequencies
!
count=count+ (focc(iband1,ik)-focc(iband2,ik))
DO iw = 1, nw
!
w = wgrid(iw)
!
jdos(iw) = jdos(iw) + (focc(iband1,ik)-focc(iband2,ik)) * &
exp(-(etrans-w)**2/intersmear**2) &
/ (intersmear * sqrt(PI))
ENDDO
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO kpt_gauss
ELSE
CALL errore('epsilon', 'invalid SMEARTYPE = '//trim(smeartype), 1)
ENDIF
!
! jdos normalizzation
!
jdos(:)=jdos(:)/count
!
! check jdos normalization
!
renorm = alpha * sum( jdos(:) )
!
! write results on data files
!
IF (ionode) THEN
WRITE(stdout,"(/,5x, 'Integration over JDOS gives: ',f15.9,' instead of 1.0d0' )") renorm
WRITE(stdout,"(/,5x, 'Writing output on file...' )")
OPEN (30, FILE='jdos.dat', FORM='FORMATTED' )
!
WRITE(30, "(2x,'# energy grid [eV] JDOS [1/eV] ')" )
!
DO iw =1, nw
!
WRITE(30,"(4f15.6)") wgrid(iw), jdos(iw)
!
ENDDO
!
CLOSE(30)
ENDIF
!
! local cleaning
!
DEALLOCATE ( jdos )
!
! collinear spin calculation
!
ELSEIF(nspin==2) THEN
!
! allocate main spectral and auxiliary quantities
!
ALLOCATE( srjdos(0:1,nw), STAT=ierr )
IF (ierr/=0) CALL errore('epsilon','allocating spin resolved jdos',abs(ierr))
!
! initialize jdos
!
srjdos(:,:)=0.0_DP
! Initialising a counter for the number of transition
srcount(:)=0.0_DP
!
! main kpt loop
!
IF (smeartype=='lorentz') THEN
DO is=0,1
! if nspin=2 the number of nks must be even (even if the calculation
! is performed at gamma point only), so nks must be always a multiple of 2
DO ik = 1 + is * int(nks/2), int(nks/2) + is * int(nks/2)
!
! Calculation of joint density of states
! 'intersmear' is the brodening parameter
!
DO iband2 = 1,nbnd
IF ( focc(iband2,ik) < 2.0d0) THEN
DO iband1 = 1,nbnd
!
IF ( focc(iband1,ik) >= 1.0d-4 ) THEN
!
! transition energy
!
etrans = ( et(iband2,ik) -et(iband1,ik) ) * RYTOEV + shift
!
IF( etrans < 1.0d-10 ) CYCLE
! loop over frequencies
!
srcount(is)=srcount(is)+ (focc(iband1,ik)-focc(iband2,ik))
DO iw = 1, nw
!
w = wgrid(iw)
!
srjdos(is,iw) = srjdos(is,iw) + intersmear * (focc(iband1,ik)-focc(iband2,ik)) &
/ ( PI * ( (etrans -w )**2 + (intersmear)**2 ) )
ENDDO
ENDIF
ENDDO
ENDIF
ENDDO
ENDDO
ENDDO
ELSEIF (smeartype=='gauss') THEN
DO is=0,1
! if nspin=2 the number of nks must be even (even if the calculation
! is performed at gamma point only), so nks must be always a multiple of 2
DO ik = 1 + is * int(nks/2), int(nks/2) + is * int(nks/2)
!
! Calculation of joint density of states
! 'intersmear' is the brodening parameter
!
DO iband2 = 1,nbnd
DO iband1 = 1,nbnd
!
IF ( focc(iband2,ik) < 2.0d0) THEN
IF ( focc(iband1,ik) >= 1.0d-4 ) THEN
!
! transition energy
!
etrans = ( et(iband2,ik) -et(iband1,ik) ) * RYTOEV + shift
!
IF( etrans < 1.0d-10 ) CYCLE
! loop over frequencies
!
srcount(is)=srcount(is)+ (focc(iband1,ik)-focc(iband2,ik))
DO iw = 1, nw
!
w = wgrid(iw)
!
srjdos(is,iw) = srjdos(is,iw) + (focc(iband1,ik)-focc(iband2,ik)) * &
exp(-(etrans-w)**2/intersmear**2) &
/ (intersmear * sqrt(PI))
ENDDO
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
ELSE
CALL errore('epsilon', 'invalid SMEARTYPE = '//trim(smeartype), 1)
ENDIF
!
! jdos normalizzation
!
DO is = 0,1
srjdos(is,:)=srjdos(is,:)/srcount(is)
ENDDO
!
! check jdos normalization
!
renormzero = alpha * sum( srjdos(0,:) )
renormuno = alpha * sum( srjdos(1,:) )
!
! write results on data files
!
IF (ionode) THEN
WRITE(stdout,"(/,5x, 'Integration over spin UP JDOS gives: ',f15.9,' instead of 1.0d0' )") renormzero
WRITE(stdout,"(/,5x, 'Integration over spin DOWN JDOS gives: ',f15.9,' instead of 1.0d0' )") renormuno
WRITE(stdout,"(/,5x, 'Writing output on file...' )")
OPEN (30, FILE='jdos.dat', FORM='FORMATTED' )
!
WRITE(30, "(2x,'# energy grid [eV] UJDOS [1/eV] DJDOS[1:eV]')" )
!
DO iw =1, nw
!
WRITE(30,"(4f15.6)") wgrid(iw), srjdos(0,iw), srjdos(1,iw)
!
ENDDO
!
CLOSE(30)
ENDIF
DEALLOCATE ( srjdos )
ENDIF
!
! local cleaning
!
CALL grid_destroy()
END SUBROUTINE jdos_calc
!-----------------------------------------------------------------------------
SUBROUTINE offdiag_calc ( intersmear,intrasmear, nw, wmax, wmin, nbndmin, nbndmax,&
shift, metalcalc, nspin )
!-----------------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE constants, ONLY : PI, RYTOEV
USE cell_base, ONLY : tpiba2, omega
USE wvfct, ONLY : nbnd, et
USE ener, ONLY : efermi => ef
USE klist, ONLY : nks, nkstot, degauss
USE grid_module, ONLY : focc, wgrid, grid_build, grid_destroy
USE io_global, ONLY : ionode, stdout
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm
USE mp, ONLY : mp_sum
!
IMPLICIT NONE
!
! input variables
!
INTEGER, INTENT(in) :: nw,nbndmin,nbndmax,nspin
REAL(DP), INTENT(in) :: wmax, wmin, intersmear,intrasmear, shift
LOGICAL, INTENT(in) :: metalcalc
!
! local variables
!
INTEGER :: ik, iband1, iband2
INTEGER :: iw, ierr, it1, it2
REAL(DP) :: etrans, const, w
!
COMPLEX(DP), ALLOCATABLE :: dipole_aux(:,:,:)
COMPLEX(DP), ALLOCATABLE :: epstot(:,:,:),dipoletot(:,:,:,:)
!
!--------------------------
! main routine body
!--------------------------
!
! perform some consistency checks, calculate occupation numbers and setup w grid
!
CALL grid_build(nw, wmax, wmin )
!
! allocate main spectral and auxiliary quantities
!
ALLOCATE( dipoletot(3,3, nbnd, nbnd), STAT=ierr )
IF (ierr/=0) CALL errore('epsilon','allocating dipoletot', abs(ierr) )
!
ALLOCATE( dipole_aux(3, nbnd, nbnd), STAT=ierr )
IF (ierr/=0) CALL errore('epsilon','allocating dipole_aux', abs(ierr) )
!
ALLOCATE(epstot( 3,3, nw),STAT=ierr )
IF (ierr/=0) CALL errore('epsilon','allocating epstot', abs(ierr))
!
! initialize response functions
!
epstot = (0.0_DP,0.0_DP)
!
! main kpt loop
!
DO ik = 1, nks
!
! For every single k-point: order k+G for
! read and distribute wavefunctions
! compute dipole matrix 3 x nbnd x nbnd parallel over g
! recover g parallelism getting the total dipole matrix
!
CALL dipole_calc( ik, dipole_aux, metalcalc, nbndmin, nbndmax)
!
DO it2 = 1, 3
DO it1 = 1, 3
dipoletot(it1,it2,:,:) = tpiba2 * dipole_aux(it1,:,:) * conjg( dipole_aux(it2,:,:) )
ENDDO
ENDDO
!
! Calculation of real and immaginary parts
! of the macroscopic dielettric function from dipole
! approximation.
! 'intersmear' is the brodening parameter
!
DO iband2 = 1,nbnd
IF ( focc(iband2,ik) < 2.0d0) THEN
DO iband1 = 1,nbnd
!
IF ( focc(iband1,ik) >= 1e-4 ) THEN
!
! transition energy
!
etrans = ( et(iband2,ik) -et(iband1,ik) ) * RYTOEV + shift
!
IF (abs(focc(iband2,ik)-focc(iband1,ik))< 1e-4) CYCLE
!
! loop over frequencies
!
DO iw = 1, nw
!
w = wgrid(iw)
!
epstot(:,:,iw) = epstot(:,:,iw) + dipoletot(:,:,iband1,iband2)*RYTOEV**3/(etrans) *&
focc(iband1,ik)/(etrans**2 - w**2 - (0,1)*intersmear*w)
ENDDO
!
ENDIF
ENDDO
ENDIF
ENDDO
!
!Intraband (only if metalcalc is true)
!
IF (metalcalc) THEN
DO iband1 = 1,nbnd
!
IF ( focc(iband1,ik) < 2.0d0) THEN
IF ( focc(iband1,ik) >= 1e-4 ) THEN
!
! loop over frequencies
!
DO iw = 1, nw
!
w = wgrid(iw)
!
epstot(:,:,iw) = epstot(:,:,iw) - dipoletot(:,:,iband1,iband1)* &
RYTOEV**2 * (exp((et(iband1,ik)-efermi)/degauss ))/ &
(( w**2 + (0,1)*intrasmear*w)*(1+exp((et(iband1,ik)-efermi)/ &
degauss))**2*degauss )
ENDDO
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
!
! recover over kpt parallelization (inter_pool)
!
CALL mp_sum( epstot, inter_pool_comm )
!
! impose the correct normalization
!
const = 64.0d0 * PI / ( omega * REAL(nkstot, DP) )
epstot(:,:,:) = epstot(:,:,:) * const
!
! add diagonal term
!
epstot(1,1,:) = 1.0_DP + epstot(1,1,:)
epstot(2,2,:) = 1.0_DP + epstot(2,2,:)
epstot(3,3,:) = 1.0_DP + epstot(3,3,:)
!
! write results on data files
!
IF (ionode) THEN
!
WRITE(stdout,"(/,5x, 'Writing output on file...' )")
!
OPEN (41, FILE='epsxx.dat', FORM='FORMATTED' )
OPEN (42, FILE='epsxy.dat', FORM='FORMATTED' )
OPEN (43, FILE='epsxz.dat', FORM='FORMATTED' )
OPEN (44, FILE='epsyx.dat', FORM='FORMATTED' )
OPEN (45, FILE='epsyy.dat', FORM='FORMATTED' )
OPEN (46, FILE='epsyz.dat', FORM='FORMATTED' )
OPEN (47, FILE='epszx.dat', FORM='FORMATTED' )
OPEN (48, FILE='epszy.dat', FORM='FORMATTED' )
OPEN (49, FILE='epszz.dat', FORM='FORMATTED' )
!
WRITE(41, "(2x,'# energy grid [eV] epsr epsi')" )
WRITE(42, "(2x,'# energy grid [eV] epsr epsi')" )
WRITE(43, "(2x,'# energy grid [eV] epsr epsi')" )
WRITE(44, "(2x,'# energy grid [eV] epsr epsi')" )
WRITE(45, "(2x,'# energy grid [eV] epsr epsi')" )
WRITE(46, "(2x,'# energy grid [eV] epsr epsi')" )
WRITE(47, "(2x,'# energy grid [eV] epsr epsi')" )
WRITE(48, "(2x,'# energy grid [eV] epsr epsi')" )
WRITE(49, "(2x,'# energy grid [eV] epsr epsi')" )
!
DO iw =1, nw
!
WRITE(41,"(4f15.6)") wgrid(iw), REAL(epstot(1,1, iw)), aimag(epstot(1,1, iw))
WRITE(42,"(4f15.6)") wgrid(iw), REAL(epstot(1,2, iw)), aimag(epstot(1,2, iw))
WRITE(43,"(4f15.6)") wgrid(iw), REAL(epstot(1,3, iw)), aimag(epstot(1,3, iw))
WRITE(44,"(4f15.6)") wgrid(iw), REAL(epstot(2,1, iw)), aimag(epstot(2,1, iw))
WRITE(45,"(4f15.6)") wgrid(iw), REAL(epstot(2,2, iw)), aimag(epstot(2,2, iw))
WRITE(46,"(4f15.6)") wgrid(iw), REAL(epstot(2,3, iw)), aimag(epstot(2,3, iw))
WRITE(47,"(4f15.6)") wgrid(iw), REAL(epstot(3,1, iw)), aimag(epstot(3,1, iw))
WRITE(48,"(4f15.6)") wgrid(iw), REAL(epstot(3,2, iw)), aimag(epstot(3,2, iw))
WRITE(49,"(4f15.6)") wgrid(iw), REAL(epstot(3,3, iw)), aimag(epstot(3,3, iw))
!
ENDDO
!
CLOSE(30)
CLOSE(40)
CLOSE(41)
CLOSE(42)
!
ENDIF
!
! local cleaning
!
CALL grid_destroy()
DEALLOCATE ( dipoletot, dipole_aux, epstot )
END SUBROUTINE offdiag_calc
!-------------------------------------------------
SUBROUTINE occ_calc ()
!-------------------------------------------------
!
USE kinds, ONLY : DP
USE klist, ONLY : nkstot, wk, degauss
USE wvfct, ONLY : nbnd, wg, et
USE ener, ONLY : ef
USE mp_global, ONLY : me_pool
!
IMPLICIT NONE
!
REAL(DP), ALLOCATABLE :: focc(:,:),foccp(:,:)
CHARACTER(25) :: filename
INTEGER :: ierr, i, ik
!
ALLOCATE ( focc( nbnd, nkstot), STAT=ierr )
IF (ierr/=0) CALL errore('grid_build','allocating focc', abs(ierr))
!
ALLOCATE ( foccp( nbnd, nkstot), STAT=ierr )
IF (ierr/=0) CALL errore('grid_build','allocating foccp', abs(ierr))
IF (me_pool==0) THEN
!
filename = 'occupations.dat'
! WRITE(filename,"(I3,'.occupation.dat')")me_pool
OPEN (unit=50, file=trim(filename))
WRITE(50,*) '#energy (Ry) occupation factor derivative'
DO ik = 1,nkstot
DO i = 1,nbnd
focc(i,ik)= wg(i, ik ) * 2.0_DP/wk( ik )
foccp(i,ik)= 2* exp((et(i,ik)-ef)/degauss)/((1+exp((et(i,ik)-ef)/degauss))**2*degauss)
WRITE(50,*)et(i,ik),focc(i,ik),foccp(i,ik)
ENDDO
ENDDO
CLOSE (50)
!
ENDIF
!
DEALLOCATE ( focc, STAT=ierr)
CALL errore('grid_destroy','deallocating grid stuff',abs(ierr))
!
DEALLOCATE ( foccp, STAT=ierr)
CALL errore('grid_destroy','deallocating grid stuff',abs(ierr))
END SUBROUTINE occ_calc
!--------------------------------------------------------------------
SUBROUTINE dipole_calc( ik, dipole_aux, metalcalc, nbndmin, nbndmax )
!------------------------------------------------------------------
!
! Evaluate overlaps (n, k-q | n', k) where k is fixed
! (it should be the first k-point in the pwscf input file),
! and q (labeled by ik) varies in the BZ
!
USE kinds, ONLY : DP
USE wvfct, ONLY : npw, nbnd, igk, g2kin, ecutwfc
USE wavefunctions_module, ONLY : evc
USE klist, ONLY : xk
USE cell_base, ONLY : tpiba2
USE gvect, ONLY : ngm, g
USE io_files, ONLY : nwordwfc, iunwfc
USE grid_module, ONLY : focc
USE mp_global, ONLY : intra_pool_comm
USE mp, ONLY : mp_sum
IMPLICIT NONE
!
! global variables
INTEGER, INTENT(in) :: ik,nbndmin,nbndmax
COMPLEX(DP), INTENT(inout) :: dipole_aux(3,nbnd,nbnd)
LOGICAL, INTENT(in) :: metalcalc
!
! local variables
INTEGER :: iband1,iband2,ig
COMPLEX(DP) :: caux
!
! Routine Body
!
CALL start_clock( 'dipole_calc' )
!
! setup k+G grids for each kpt
!
CALL gk_sort (xk (1, ik), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
!
! read wfc for the given kpt
!
CALL davcio (evc, nwordwfc, iunwfc, ik, - 1)
!
! compute matrix elements
!
dipole_aux(:,:,:) = (0.0_DP,0.0_DP)
!
DO iband2 = nbndmin,nbndmax
IF ( focc(iband2,ik) < 2.0d0) THEN
DO iband1 = nbndmin,nbndmax
!
IF ( iband1==iband2 ) CYCLE
IF ( focc(iband1,ik) >= 1e-4 ) THEN
!
DO ig=1,npw
!
caux= conjg(evc(ig,iband1))*evc(ig,iband2)
!
dipole_aux(:,iband1,iband2) = dipole_aux(:,iband1,iband2) + &
( g(:,igk(ig)) ) * caux
!
ENDDO
ENDIF
!
ENDDO
ENDIF
ENDDO
!
! The diagonal terms are taken into account only if the system is treated like a metal, not
! in the intraband therm. Because of this we can recalculate the diagonal component of the dipole
! tensor directly as we need it for the intraband therm, without interference with interband one.
!
IF (metalcalc) THEN
!
DO iband1 = nbndmin,nbndmax
DO ig=1,npw
!
caux= conjg(evc(ig,iband1))*evc(ig,iband1)
!
dipole_aux(:,iband1,iband1) = dipole_aux(:,iband1,iband1) + &
( g(:,igk(ig))+ xk(:,ik) ) * caux
!
ENDDO
ENDDO
!
ENDIF
!
! recover over G parallelization (intra_pool)
!
CALL mp_sum( dipole_aux, intra_pool_comm )
!
CALL stop_clock( 'dipole_calc' )
!
END SUBROUTINE dipole_calc
!
!----------------------------------------------------------------
SUBROUTINE refold ( gmap, nshell, g0vec )
!----------------------------------------------------------------
!
! the map is defined as follows
! g(:,gmap(ig,ig0)) = g(:,ig) - g0vec(:,ig0)
!
! at the exit, the folding vectors are in cartesian coordinates
!
!----------------------------------------------------------------
!
USE parameters
USE gvect !, ONLY: g, ngm
USE cell_base, ONLY: bg, at
USE io_global, ONLY: ionode
implicit none
!
integer :: notfound, g2(3), count
integer :: ig0, ig, igp, ig1, ig2, ig3
integer :: nshell, ng0vec, npw
real(DP) :: g1(3)
real(DP) :: g0vec(3,(nshell*2+1)**3)
integer :: gmap(ngm,(nshell*2+1)**3)
real(DP) :: gtmp(3,ngm)
!
! the 3^3 possible translations (in very anisotropic materials
! this could go to 5^3 or even more - see EPW code)
!
! crystal coordinates
!
ng0vec = 0
do ig1 = -nshell,nshell
do ig2 = -nshell,nshell
do ig3 = -nshell,nshell
ng0vec = ng0vec + 1
!print*, ng0vec
g0vec(1,ng0vec) = ig1
g0vec(2,ng0vec) = ig2
g0vec(3,ng0vec) = ig3
enddo
enddo
enddo
!
! bring all the G-vectors in crystal coordinates
! (in real codes we use the array of the miller indices mill_ )
!
! gtmp(:,:) = g(:,:)
! call cryst_to_cart (ngm, gtmp, at, -1)
call cryst_to_cart (ngm, g, at, -1)
!
count = 0
notfound = 0
do ig = 1, ngm
!
do ig0 = 1, ng0vec
!
gmap (ig,ig0) = 0
count = count + 1
!
g2 = nint ( g(:,ig) - g0vec(:,ig0) )
!
do igp = 1, ngm
! if ( ( g2(1) .eq. nint(gtmp(1,igp)) ) .and. &
! ( g2(2) .eq. nint(gtmp(2,igp)) ) .and. &
! ( g2(3) .eq. nint(gtmp(3,igp)) ) ) &
if ( ( g2(1) .eq. nint(g(1,igp)) ) .and. &
( g2(2) .eq. nint(g(2,igp)) ) .and. &
( g2(3) .eq. nint(g(3,igp)) ) ) &
gmap (ig,ig0) = igp
! if(ionode)print*,' '
! if(ionode)print*,' --- ', ig, ig0 ,igp, gmap(ig,ig0),' --- '
! if(ionode)print*,g2(:)
! if(ionode)print*,g(:,igp)
! if(ionode)print*,g(:,gmap(ig,ig0))
enddo
if (gmap (ig,ig0).eq.0) notfound = notfound + 1
!
enddo
!
enddo
write(6,'(4x,"refold: notfound = ",i6," out of ",i6)') notfound,count
!if(ionode) print*, gmap
!
! back to cartesian
!
call cryst_to_cart (ngm, g, bg, 1)
call cryst_to_cart (ng0vec, g0vec, bg, 1)
!
END SUBROUTINE refold
!
!
!----------------------------------------------------------------
SUBROUTINE refold2 ( gmap, nshell, g0vec )
!----------------------------------------------------------------
!
! the map is defined as follows
! g(:,gmap(ig,ig0)) = g(:,ig) - g0vec(:,ig0)
!
! at the exit, the folding vectors are in cartesian coordinates
!
!----------------------------------------------------------------
!
USE parameters
USE gvect !, ONLY: g, ngm
USE cell_base, ONLY: bg, at
USE io_global, ONLY: ionode
implicit none
!
integer :: notfound, g2(3), count
integer :: ig0, ig, igp, ig1, ig2, ig3
integer :: nshell, ng0vec, npw
! parameter :: ng0vec
real(DP) :: g1(3)
real(DP) :: g0vec(3,(nshell*2+1)**3)
integer :: gmap(ngm,(nshell*2+1)**3)
!
! the 3^3 possible translations (in very anisotropic materials
! this could go to 5^3 or even more - see EPW code)
!
! crystal coordinates
!
ng0vec = 0
do ig1 = -nshell,nshell
do ig2 = -nshell,nshell
do ig3 = -nshell,nshell
ng0vec = ng0vec + 1
!print*, ng0vec
g0vec(1,ng0vec) = ig1
g0vec(2,ng0vec) = ig2
g0vec(3,ng0vec) = ig3
enddo
enddo
enddo
!
! bring all the G-vectors in crystal coordinates
! (in real codes we use the array of the miller indices mill_ )
!
! call cryst_to_cart (ngm, g, at, -1)
!
count = 0
notfound = 0
do ig = 1, ngm
! if(ionode)print*,g(:,ig)
!
do ig0 = 1, ng0vec
!
gmap (ig,ig0) = 0
count = count + 1
!
g2 = nint( g(:,ig) - g0vec(:,ig0) )
!
do igp = 1, ngm
if ( ( g2(1) .eq. nint(g(1,igp)) ) .and. &
( g2(2) .eq. nint(g(2,igp)) ) .and. &
( g2(3) .eq. nint(g(3,igp)) ) ) then
gmap (ig,ig0) = igp
! if(ionode)print*,' '
! if(ionode)print*,' --- ', ig, ig0 ,igp, gmap(ig,ig0),' --- '
! if(ionode)print*,g2(:)
! if(ionode)print*,g(:,igp)
! if(ionode)print*,g(:,gmap(ig,ig0))
endif
enddo
if (gmap (ig,ig0).eq.0) notfound = notfound + 1
!
enddo
!
enddo
write(6,'(4x,"refold: notfound = ",i6," out of ",i6)') notfound,count
!if(ionode) print*, gmap
!
! back to cartesian
!
! call cryst_to_cart (ngm, g, bg, 1)
! call cryst_to_cart (ng0vec, g0vec, bg, 1)
!
END SUBROUTINE refold2
!
SUBROUTINE get_vc( rho, rho_core, rhog_core, etxc, vtxc, v )
!----------------------------------------------------------------------------
!
! ... Exchange-Correlation potential Vxc(r) from n(r)
!
USE kinds, ONLY : DP
USE constants, ONLY : e2, eps8
USE io_global, ONLY : stdout, ionode
USE fft_base, ONLY : dfftp
USE gvect, ONLY : ngm
USE lsda_mod, ONLY : nspin
USE cell_base, ONLY : omega
USE spin_orb, ONLY : domag
USE funct, ONLY : xc, xc_spin
USE scf, ONLY : scf_type
USE mp_global, ONLY : intra_pool_comm, intra_bgrp_comm, mpime
USE mp, ONLY : mp_sum
!
IMPLICIT NONE
!
TYPE (scf_type), INTENT(IN) :: rho
REAL(DP), INTENT(IN) :: rho_core(dfftp%nnr)
! the core charge
COMPLEX(DP), INTENT(IN) :: rhog_core(ngm)
! input: the core charge in reciprocal space
REAL(DP), INTENT(OUT) :: v(dfftp%nnr,nspin), vtxc, etxc
! V_xc potential
! integral V_xc * rho
! E_xc energy
!
! ... local variables
!
REAL(DP) :: rhox, arhox, zeta, amag, vs, ex, ec, vx(2), vc(2), rhoneg(2)
! the total charge in each point
! the absolute value of the charge
! the absolute value of the charge
! local exchange energy
! local correlation energy
! local exchange potential
! local correlation potential
INTEGER :: ir, ipol
! counter on mesh points
! counter on nspin
!
REAL(DP), PARAMETER :: vanishing_charge = 1.D-10, &
vanishing_mag = 1.D-20
!
!
CALL start_clock( 'v_xc' )
!
etxc = 0.D0
vtxc = 0.D0
v(:,:) = 0.D0
rhoneg = 0.D0
!
IF ( nspin == 1 .OR. ( nspin == 4 .AND. .NOT. domag ) ) THEN
!
! ... spin-unpolarized case
!
!$omp parallel do private( rhox, arhox, ex, ec, vx, vc ), &
!$omp reduction(+:etxc,vtxc), reduction(-:rhoneg)
DO ir = 1, dfftp%nnr
!
rhox = rho%of_r(ir,1) + rho_core(ir)
!
arhox = ABS( rhox )
!
IF ( arhox > vanishing_charge ) THEN
!
CALL xc( arhox, ex, ec, vx(1), vc(1) )
!
v(ir,1) = e2*( vc(1) )
!
etxc = etxc + e2*( ex + ec ) * rhox
!
vtxc = vtxc + v(ir,1) * rho%of_r(ir,1)
!
ENDIF
!
IF ( rho%of_r(ir,1) < 0.D0 ) rhoneg(1) = rhoneg(1) - rho%of_r(ir,1)
!
END DO
!$omp end parallel do
!
ELSE IF ( nspin == 2 ) THEN
! IF(ionode) THEN
! PRINT*, 'This is supposed to work only with ispin = 1 '
! ENDIF
! stop
! ENDIF
!
! ... spin-polarized case
!
!$omp parallel do private( rhox, arhox, zeta, ex, ec, vx, vc ), &
!$omp reduction(+:etxc,vtxc), reduction(-:rhoneg)
DO ir = 1, dfftp%nnr
!
rhox = rho%of_r(ir,1) + rho%of_r(ir,2) + rho_core(ir)
!
arhox = ABS( rhox )
!
IF ( arhox > vanishing_charge ) THEN
!
zeta = ( rho%of_r(ir,1) - rho%of_r(ir,2) ) / arhox
!
IF ( ABS( zeta ) > 1.D0 ) zeta = SIGN( 1.D0, zeta )
!
IF ( rho%of_r(ir,1) < 0.D0 ) rhoneg(1) = rhoneg(1) - rho%of_r(ir,1)
IF ( rho%of_r(ir,2) < 0.D0 ) rhoneg(2) = rhoneg(2) - rho%of_r(ir,2)
!
CALL xc_spin( arhox, zeta, ex, ec, vx(1), vx(2), vc(1), vc(2) )
!
v(ir,:) = e2*( vx(:) + vc(:) )
!
etxc = etxc + e2*( ex + ec ) * rhox
!
vtxc = vtxc + v(ir,1) * rho%of_r(ir,1) + v(ir,2) * rho%of_r(ir,2)
!
END IF
!
END DO
!$omp end parallel do
!
ELSE IF ( nspin == 4 ) THEN
!
! ... noncolinear case
!
DO ir = 1,dfftp%nnr
!
amag = SQRT( rho%of_r(ir,2)**2 + rho%of_r(ir,3)**2 + rho%of_r(ir,4)**2 )
!
rhox = rho%of_r(ir,1) + rho_core(ir)
!
IF ( rho%of_r(ir,1) < 0.D0 ) rhoneg(1) = rhoneg(1) - rho%of_r(ir,1)
!
arhox = ABS( rhox )
!
IF ( arhox > vanishing_charge ) THEN
!
zeta = amag / arhox
!
IF ( ABS( zeta ) > 1.D0 ) THEN
!
rhoneg(2) = rhoneg(2) + 1.D0 / omega
!
zeta = SIGN( 1.D0, zeta )
!
END IF
!
CALL xc_spin( arhox, zeta, ex, ec, vx(1), vx(2), vc(1), vc(2) )
!
vs = 0.5D0*( vx(1) + vc(1) - vx(2) - vc(2) )
!
v(ir,1) = e2*( 0.5D0*( vx(1) + vc(1) + vx(2) + vc(2 ) ) )
!
IF ( amag > vanishing_mag ) THEN
!
DO ipol = 2, 4
!
v(ir,ipol) = e2 * vs * rho%of_r(ir,ipol) / amag
!
vtxc = vtxc + v(ir,ipol) * rho%of_r(ir,ipol)
!
END DO
!
END IF
!
etxc = etxc + e2*( ex + ec ) * rhox
vtxc = vtxc + v(ir,1) * rho%of_r(ir,1)
!
END IF
!
END DO
!
END IF
!
CALL mp_sum( rhoneg , intra_bgrp_comm )
!
rhoneg(:) = rhoneg(:) * omega / ( dfftp%nr1*dfftp%nr2*dfftp%nr3 )
!
IF ( rhoneg(1) > eps8 .OR. rhoneg(2) > eps8 ) &
WRITE( stdout,'(/,5X,"negative rho (up, down): ",2E10.3)') rhoneg
!
! ... energy terms, local-density contribution
!
vtxc = omega * vtxc / ( dfftp%nr1*dfftp%nr2*dfftp%nr3 )
etxc = omega * etxc / ( dfftp%nr1*dfftp%nr2*dfftp%nr3 )
!
! ... add gradient corrections (if any)
!
! CALL gradcorr( rho%of_r, rho%of_g, rho_core, rhog_core, etxc, vtxc, v )
!
! ... add non local corrections (if any)
!
! CALL nonloccorr(rho%of_r, rho_core, etxc, vtxc, v)
!
CALL mp_sum( vtxc , intra_bgrp_comm )
CALL mp_sum( etxc , intra_bgrp_comm )
!
CALL stop_clock( 'v_xc' )
!
RETURN
!
END SUBROUTINE get_vc
SUBROUTINE GN_ppa (nw,eta, wgrid, ieps, wtilde, wplasma , sysd)
USE kinds, ONLY : DP
USE io_global, ONLY : stdout, ionode
integer nw ! frequency index for the fitting of rht PPA parameter
integer sysd ! dimensionality of the system
real(DP) wgrid (nw) ! (imaginary) frequency for the PPA parameter
REAL(DP) :: ieps(3,nw)
real(DP) wtilde, wplasma ! PPA parameters
real(DP) eta
!---- internal variables ------
real(DP) eps0 ,eps1
integer iw
!use Godby Needs plasmon-pole model
!for 3D
if (sysd == 3)then
iw = nw/2
eps0 =1.d0/(sum(ieps(:,1))/3.d0)
eps1 =1.d0/(sum(ieps(:,iw))/3.d0)
elseif(sysd ==2 )then
iw = nw/2
eps0 =1.d0/((ieps(1,1) +ieps(1,1) )/2.d0)
eps1 =1.d0/((ieps(1,iw)+ieps(2,iw))/2.d0)
endif
wtilde = (( eps0 - 1.d0 )/(eps0-eps1) - 1.d0) * wgrid(iw)**2
if (wtilde.gt.0)then
wtilde = sqrt(wtilde)
else
if(ionode)print*, 'wtilde is imaginary! Stop!'
stop
endif
wplasma = (1.d0 - eps0)*wtilde**2
if (wplasma.gt.0)then
wplasma = sqrt(wplasma)
else
if(ionode)print*, 'wplasma is imaginary! Stop!'
stop
endif
IF(ionode) THEN
PRINT *
PRINT *, ' Godby-Needs Plasmon-pole parameters : '
PRINT *, ' wtilde = ',wtilde
PRINT *, ' wplasma = ',wplasma
PRINT *
OPEN (30, FILE='epsppa_GN.dat')
DO jw = 1, nw
WRITE(30,*) wgrid(jw), &
real(1.0d0/(1.0d0+ wplasma**2/(wgrid(jw)**2-(wtilde-(0.0d0,1.0d0)*eta/100)**2))) , &
aimag(1.0d0/(1.0d0+ wplasma**2/(wgrid(jw)**2-(wtilde-(0.0d0,1.0d0)*eta/100)**2))) !, &
!epsr(1,jw)
ENDDO
CLOSE(30)
OPEN (30, FILE='epsppa_GN_inv.dat')
DO jw = 1, nw
WRITE(30,*) wgrid(jw), &
real((1.0d0+ wplasma**2/(wgrid(jw)**2-(wtilde-(0.0d0,1.0d0)*eta/100)**2))) , &
aimag((1.0d0+ wplasma**2/(wgrid(jw)**2-(wtilde-(0.0d0,1.0d0)*eta/100)**2))) !, &
!epsr(1,jw)
ENDDO
CLOSE(30)
ENDIF
END SUBROUTINE GN_ppa
SUBROUTINE print_epsilon (npoles, par, wgrid, nw )
USE io_global, ONLY : stdout, ionode
USE kinds, ONLY : DP
integer npoles
real(DP) par(npoles*2)
integer nw
real(DP) wgrid(nw)
integer iw, ipole
real eta
complex(DP) epsinv
eta = 0.01
IF(ionode)THEN
PRINT*, ' I recreating epsilon, and printing it to file'
PRINT*, par
OPEN (30, FILE='epsilon_new.dat')
DO iw = 1, nw, 1
epsinv = (1.d0 , 0.d0)
DO ipole = 1, npoles
epsinv = epsinv + par(ipole*2)**2 / ( wgrid(iw)**2-(par(ipole*2-1) - (0.0d0,1.0d0)*eta)**2 )
ENDDO
WRITE(30,*) wgrid(iw), real (1.d0 / epsinv )
!real(1.0d0/(1.0d0+ wplasma0**2/(wgrid(iw)**2-(wtilde-(0.0d0,1.0d0)*eta)**2))), &
ENDDO
CLOSE(30)
ENDIF
END SUBROUTINE print_epsilon