Commit ff55a2b3 authored by Marco Govoni's avatar Marco Govoni
Browse files

Upgraded class bz_grid, propagated changes to all code except solve_qp and solve_gfreq.

parent d8f7441e
......@@ -65,9 +65,9 @@ SUBROUTINE store_sqvc(sqvc_tmp,numg,singularity_removal_mode,div,printout_div)
nq(2)=1._DP
nq(3)=1._DP
ELSE
nq(1) = REAL( q_grid%np1, KIND=DP )
nq(2) = REAL( q_grid%np2, KIND=DP )
nq(3) = REAL( q_grid%np3, KIND=DP )
nq(1) = REAL( q_grid%ngrid(1), KIND=DP )
nq(2) = REAL( q_grid%ngrid(2), KIND=DP )
nq(3) = REAL( q_grid%ngrid(3), KIND=DP )
ENDIF
!
! ======
......@@ -322,11 +322,11 @@ SUBROUTINE store_sqvc_q(sqvc_tmp,numg,singularity_removal_mode,iq,l_use_igq)
!
CALL start_clock( 'storesqvcq' )
!
nq(1) = REAL( q_grid%np1, KIND=DP )
nq(2) = REAL( q_grid%np2, KIND=DP )
nq(3) = REAL( q_grid%np3, KIND=DP )
nq(1) = REAL( q_grid%ngrid(1), KIND=DP )
nq(2) = REAL( q_grid%ngrid(2), KIND=DP )
nq(3) = REAL( q_grid%ngrid(3), KIND=DP )
!
q(:) = q_grid%xp_cart(:,iq)
q(:) = q_grid%p_cart(:,iq)
!
! =======
! BODY
......
......@@ -359,7 +359,7 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
USE io_push, ONLY : io_push_title
USE mp_world, ONLY : mpime,world_comm
USE class_bz_grid, ONLY : bz_grid
USE types_bz_grid, ONLY : kmq_grid,q_grid
USE types_bz_grid, ONLY : q_grid, compute_phase
!
IMPLICIT NONE
!
......@@ -373,7 +373,7 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
!
! Workspace
!
INTEGER :: ipert, ig, ir, ibnd, iks, ikqs
INTEGER :: ipert, ig, ir, ibnd, iks, ikqs, ik
! Counter on perturbations
! Counter on plane waves
! Counter on real space grids
......@@ -386,6 +386,7 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
! Current k-q point
INTEGER :: npwkq
!
REAL(DP) :: kmq(3), g0(3)
REAL(DP) :: anorm, prod
REAL(DP), ALLOCATABLE :: eprec(:)
! Preconditioning matrix
......@@ -426,16 +427,21 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
ALLOCATE( phase(dffts%nnr) )
!
DO iks = 1, nks ! KPOINT-SPIN LOOP
!
ik = k_grid%ip(iks)
!
current_k = iks
!
ikqs = kmq_grid%index_kq(iks,iq)
k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
ikqs = k_grid%find( kmq, 'cart' )
!ikqs = kmq_grid%index_kq(iks,iq)
!
! computes the phase needed to bring the wavefunction at k-q
! to the equivalent [k-q] point in the first BZ
!
CALL kmq_grid%get_phase(iks,iq)
phase = kmq_grid%phase
CALL compute_phase( g0, 'cart', phase )
!CALL kmq_grid%get_phase(iks,iq)
!phase = kmq_grid%phase
!
IF ( lsda ) current_spin = isk(iks)
CALL g2_kin(iks)
......@@ -453,7 +459,7 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
!
! ... More stuff needed by the hamiltonian: nonlocal projectors
!
IF ( nkb > 0 ) CALL init_us_2( ngk(iks), igk_k(1,iks), xk(1,iks), vkb )
IF ( nkb > 0 ) CALL init_us_2( ngk(iks), igk_k(1,iks), k_grid%p_cart(1,ik), vkb )
!
! Read wavefuctions at k in G space, for all bands,
! and store them in evc
......@@ -632,7 +638,7 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
!
ENDDO ! K-POINT and SPIN
!
IF ( q_grid%l_gammap(iq) ) THEN
IF ( q_grid%l_pIsGamma(iq) ) THEN
IF ( gstart == 2 ) dng(1,1:m) = CMPLX( 0._DP, 0._DP, KIND=DP )
ENDIF
!
......@@ -646,3 +652,4 @@ SUBROUTINE dfpt_q (m,dvg,dng,tr2,iq)
CALL stop_bar_type( barra, 'dfpt_q' )
!
END SUBROUTINE
......@@ -557,8 +557,8 @@ MODULE wfreq_restart
!
TYPE(bks_type),INTENT(IN) :: bks
INTEGER, INTENT(IN) :: npg, npl
COMPLEX(DP), INTENT(IN) :: dmat(npg,npl,ifr%nloc,q_grid%nps)
COMPLEX(DP), INTENT(IN) :: zmat(npg,npl,rfr%nloc,q_grid%nps)
COMPLEX(DP), INTENT(IN) :: dmat(npg,npl,ifr%nloc,q_grid%np)
COMPLEX(DP), INTENT(IN) :: zmat(npg,npl,rfr%nloc,q_grid%np)
!
! Workspace
!
......@@ -585,9 +585,9 @@ MODULE wfreq_restart
!
! DMAT
!
ALLOCATE( tmp_dmat( npg, npl, ifr%nglob, q_grid%nps ))
ALLOCATE( tmp_dmat( npg, npl, ifr%nglob, q_grid%np ))
tmp_dmat = 0._DP
DO iq = 1, q_grid%nps
DO iq = 1, q_grid%np
DO ip = 1, ifr%nloc
ip_glob = ifr%l2g(ip)
tmp_dmat(:,:,ip_glob,iq) = dmat(:,:,ip,iq)
......@@ -600,15 +600,15 @@ MODULE wfreq_restart
fname=TRIM( dirname ) // '/' // TRIM(my_label)
iunit = 2000 + my_image_id
lproc = (me_bgrp == 0)
CALL serial_data_write(lproc,iunit,fname,tmp_dmat,npg,npl,ifr%nglob,q_grid%nps)
CALL serial_data_write(lproc,iunit,fname,tmp_dmat,npg,npl,ifr%nglob,q_grid%np)
!
DEALLOCATE(tmp_dmat)
!
! ZMAT
!
ALLOCATE( tmp_zmat( npg, npl, rfr%nglob, q_grid%nps ))
ALLOCATE( tmp_zmat( npg, npl, rfr%nglob, q_grid%np ))
tmp_zmat = 0._DP
DO iq = 1, q_grid%nps
DO iq = 1, q_grid%np
DO ip = 1, rfr%nloc
ip_glob = rfr%l2g(ip)
tmp_zmat(:,:,ip_glob,iq) = zmat(:,:,ip,iq)
......@@ -621,7 +621,7 @@ MODULE wfreq_restart
fname=TRIM( dirname ) // '/' // TRIM(my_label)
iunit = 2000 + my_image_id
lproc = (me_bgrp == 0)
CALL serial_data_write(lproc,iunit,fname,tmp_zmat,npg,npl,rfr%nglob,q_grid%nps)
CALL serial_data_write(lproc,iunit,fname,tmp_zmat,npg,npl,rfr%nglob,q_grid%np)
!
DEALLOCATE(tmp_zmat)
!
......@@ -869,8 +869,8 @@ MODULE wfreq_restart
!
TYPE(bks_type), INTENT(OUT) :: bks
INTEGER, INTENT(IN) :: npg, npl
COMPLEX(DP), INTENT(OUT) :: dmat(npg,npl,ifr%nloc,q_grid%nps)
COMPLEX(DP), INTENT(OUT) :: zmat(npg,npl,rfr%nloc,q_grid%nps)
COMPLEX(DP), INTENT(OUT) :: dmat(npg,npl,ifr%nloc,q_grid%np)
COMPLEX(DP), INTENT(OUT) :: zmat(npg,npl,rfr%nloc,q_grid%np)
!
! Workspace
!
......@@ -906,17 +906,17 @@ MODULE wfreq_restart
!
! READ
!
ALLOCATE( tmp_dmat( npg, npl, ifr%nglob, q_grid%nps ))
ALLOCATE( tmp_dmat( npg, npl, ifr%nglob, q_grid%np ))
!
WRITE(my_label,'("dmat_iq",i5.5,"_iks",i5.5,"_iv",i5.5,"_I",i6.6)') bks%lastdone_q,bks%lastdone_ks,bks%lastdone_band,&
& my_image_id
fname=TRIM( dirname ) // '/' // TRIM(my_label)
lproc = (me_bgrp==0)
iunit = 2000 + my_image_id
CALL serial_data_read(lproc,iunit,fname,tmp_dmat,npg,npl,ifr%nglob,q_grid%nps)
CALL serial_data_read(lproc,iunit,fname,tmp_dmat,npg,npl,ifr%nglob,q_grid%np)
!
CALL mp_bcast( tmp_dmat, 0, intra_bgrp_comm)
DO iq = 1, q_grid%nps
DO iq = 1, q_grid%np
DO ip = 1, ifr%nloc
ip_glob = ifr%l2g(ip)
dmat(:,:,ip,iq) = tmp_dmat(:,:,ip_glob,iq)
......@@ -924,17 +924,17 @@ MODULE wfreq_restart
ENDDO
DEALLOCATE(tmp_dmat)
!
ALLOCATE( tmp_zmat( npg, npl, rfr%nglob, q_grid%nps ))
ALLOCATE( tmp_zmat( npg, npl, rfr%nglob, q_grid%np ))
!
WRITE(my_label,'("zmat_iq",i5.5,"_iks",i5.5,"_iv",i5.5,"_I",i6.6)') bks%lastdone_q,bks%lastdone_ks,bks%lastdone_band,&
& my_image_id
fname=TRIM( dirname ) // '/' // TRIM(my_label)
lproc = (me_bgrp==0)
iunit = 2000 + my_image_id
CALL serial_data_read(lproc,iunit,fname,tmp_zmat,npg,npl,rfr%nglob,q_grid%nps)
CALL serial_data_read(lproc,iunit,fname,tmp_zmat,npg,npl,rfr%nglob,q_grid%np)
!
CALL mp_bcast( tmp_zmat, 0, intra_bgrp_comm)
DO iq = 1, q_grid%nps
DO iq = 1, q_grid%np
DO ip = 1, rfr%nloc
ip_glob = rfr%l2g(ip)
zmat(:,:,ip,iq) = tmp_zmat(:,:,ip_glob,iq)
......
This diff is collapsed.
......@@ -60,24 +60,27 @@ SUBROUTINE do_setup
!
CALL set_dirs()
!
! INIT K, Q GRIDS
! INIT K GRID
!
k_grid = bz_grid()
CALL k_grid%init('K')
!
IF ( .NOT. gamma_only ) THEN
!
! INIT Q GRID
!
! initialize q-point grid
q_grid = bz_grid()
CALL q_grid%init('Q')
IF (q_grid%np1 /= k_grid%np1 .OR. q_grid%np2 /= k_grid%np2 .OR. q_grid%np3 /= k_grid%np3) THEN
IF ( ANY ( q_grid%ngrid(:) - k_grid%ngrid(:) /= 0 ) ) THEN
CALL errore( 'do_setup','q-point grid must be the same as k-point grid ',1)
ENDIF
! initialize (k-q) grid
kmq_grid = bz_grid()
CALL kmq_grid%init_kq( k_grid, q_grid, -1 )
! initialize (k-q) grid
kpq_grid = bz_grid()
CALL kpq_grid%init_kq( k_grid, q_grid, +1 )
!kmq_grid = bz_grid()
!CALL kmq_grid%init_kq( k_grid, q_grid, -1 )
!! initialize (k-q) grid
!kpq_grid = bz_grid()
!CALL kpq_grid%init_kq( k_grid, q_grid, +1 )
ENDIF
!
IF( mpime == root ) THEN
......@@ -207,72 +210,45 @@ SUBROUTINE do_setup
CALL json%add('system.cell.alat',alat)
ENDIF
!
! WRITE( stdout, '(5x,"number of ks points=",i6)') nkstot
! IF( mpime == root ) CALL json%add('system.kpt.nkstot',nkstot)
! WRITE( stdout, '(23x,"cart. coord. in units 2pi/alat")')
! DO ik = 1, nkstot
! WRITE( cik, '(i6)') ik
! WRITE( stdout, '(8x,"k(",i5,") = (",3f12.7,"), wk =",f12.7)') ik, &
! (xk (ipol, ik) , ipol = 1, 3) , wk (ik)
! IF( mpime == root ) CALL json%add('system.kpt.k('//TRIM(ADJUSTL(cik))//').cartcoord:tpiba',xk(1:3))
! IF( mpime == root ) CALL json%add('system.kpt.k('//TRIM(ADJUSTL(cik))//').weight',wk(ik))
! ENDDO
! WRITE( stdout, '(/23x,"cryst. coord.")')
! DO ik = 1, nkstot
! WRITE( cik, '(i6)') ik
! DO ipol = 1, 3
! xkg(ipol) = at(1,ipol)*xk(1,ik) + at(2,ipol)*xk(2,ik) + &
! at(3,ipol)*xk(3,ik)
! ! xkg are the component in the crystal RL basis
! ENDDO
! WRITE( stdout, '(8x,"k(",i5,") = (",3f12.7,"), wk =",f12.7)') &
! ik, (xkg (ipol) , ipol = 1, 3) , wk (ik)
! IF( mpime == root ) CALL json%add('system.kpt.k('//TRIM(ADJUSTL(cik))//').crystcoord',xkg(1:3))
! ENDDO
WRITE( stdout, '(5x,"number of ks points = ",i6)') k_grid%nps
IF( mpime == root ) CALL json%add('system.kpt.nkstot',k_grid%nps)
WRITE( stdout, '(23x,"cart. coord. in units 2pi/alat")')
DO ik = 1, k_grid%nps
DO iks = 1, k_grid%nps
ik = k_grid%ip(iks)
WRITE( cik, '(i6)') ik
WRITE( stdout, '(8x,"k(",i5,") = (",3f12.7,"), wk =",f12.7)') ik, &
(k_grid%xp_cart(ipol,ik) , ipol = 1, 3) , k_grid%wp(ik)
IF( mpime == root ) CALL json%add('system.kpt.k('//TRIM(ADJUSTL(cik))//').cartcoord:tpiba',k_grid%xp_cart(1:3,ik))
IF( mpime == root ) CALL json%add('system.kpt.k('//TRIM(ADJUSTL(cik))//').weight',k_grid%wp(ik))
WRITE( stdout, '(8x,"k(",i5,") = (",3f12.7,"), wk =",f12.7)') iks, &
(k_grid%p_cart(ipol,ik) , ipol = 1, 3) , k_grid%weight(iks)
IF( mpime == root ) CALL json%add('system.kpt.k('//TRIM(ADJUSTL(cik))//').cartcoord:tpiba',k_grid%p_cart(1:3,ik))
IF( mpime == root ) CALL json%add('system.kpt.k('//TRIM(ADJUSTL(cik))//').weight',k_grid%weight(iks))
ENDDO
WRITE( stdout, '(/23x,"cryst. coord.")')
DO ik = 1, k_grid%nps
ik = k_grid%ip(iks)
WRITE( cik, '(i6)') ik
WRITE( stdout, '(8x,"k(",i5,") = (",3f12.7,"), wk =",f12.7)') &
ik, (k_grid%xp_cryst(ipol,ik) , ipol = 1, 3) , k_grid%wp(ik)
IF( mpime == root ) CALL json%add('system.kpt.k('//TRIM(ADJUSTL(cik))//').crystcoord',k_grid%xp_cryst(1:3,ik))
ik, (k_grid%p_cryst(ipol,ik) , ipol = 1, 3) , k_grid%weight(iks)
IF( mpime == root ) CALL json%add('system.kpt.k('//TRIM(ADJUSTL(cik))//').crystcoord',k_grid%p_cryst(1:3,ik))
ENDDO
!
! q-point grid
!
IF (.NOT. gamma_only ) THEN
WRITE( stdout, * )
WRITE( stdout, '(5x,"number of q points = ",i6)') q_grid%nps
IF( mpime == root ) CALL json%add('system.qpt.nqtot',q_grid%nps)
WRITE( stdout, '(5x,"number of q points = ",i6)') q_grid%np
IF( mpime == root ) CALL json%add('system.qpt.nqtot',q_grid%np)
WRITE( stdout, '(23x,"cart. coord. in units 2pi/alat")')
DO iq = 1, q_grid%nps
DO iq = 1, q_grid%np
WRITE( ciq, '(i6)') iq
WRITE( stdout, '(8x,"q(",i5,") = (",3f12.7,")")') iq, &
(q_grid%xp_cart(ipol, iq) , ipol = 1, 3)
IF( mpime == root ) CALL json%add('system.qpt.q('//TRIM(ADJUSTL(ciq))//').cartcoord:tpiba',q_grid%xp_cart(1:3,iq))
(q_grid%p_cart(ipol, iq) , ipol = 1, 3)
IF( mpime == root ) CALL json%add('system.qpt.q('//TRIM(ADJUSTL(ciq))//').cartcoord:tpiba',q_grid%p_cart(1:3,iq))
ENDDO
WRITE( stdout, '(/23x,"cryst. coord.")')
DO iq = 1, q_grid%nps
DO iq = 1, q_grid%np
WRITE( ciq, '(i6)') iq
WRITE( stdout, '(8x,"q(",i5,") = (",3f12.7,")")') &
iq, (q_grid%xp_cryst(ipol,iq) , ipol = 1, 3)
IF( mpime == root ) CALL json%add('system.qpt.q('//TRIM(ADJUSTL(ciq))//').crystcoord',q_grid%xp_cryst(1:3,iq))
ENDDO
WRITE(stdout, '(/5x,"setup uniform grid of", i3, " q-points centered on each k-point")') kmq_grid%nps
WRITE( stdout, '(5x,"(k-q)-points:")')
WRITE( stdout, '(23x,"cart. coord. in units 2pi/alat")' )
DO ik = 1, kmq_grid%nps
WRITE( stdout, '(8x,"k(",i5,") = (",3f12.7,")")') ik, &
(kmq_grid%xp_cart(ipol, ik) , ipol = 1, 3)
iq, (q_grid%p_cryst(ipol,iq) , ipol = 1, 3)
IF( mpime == root ) CALL json%add('system.qpt.q('//TRIM(ADJUSTL(ciq))//').crystcoord',q_grid%p_cryst(1:3,iq))
ENDDO
ENDIF
!
......
......@@ -42,7 +42,6 @@ SUBROUTINE read_pwout() ! to be sync'd with PW/src/read_file.f90
USE gvecw, ONLY : ecutwfc
USE fft_base, ONLY : dfftp
USE fft_base, ONLY : dffts
! USE wvfct, ONLY : npwx
USE control_flags, ONLY : gamma_only
!
IMPLICIT NONE
......
......@@ -57,16 +57,16 @@ SUBROUTINE set_npwq()
fftdriver = 'Wave'
! 'Dense' grid not yet implemented
!
npwqx = n_plane_waves( gcutw, q_grid%nps, q_grid%xp_cart, g, ngm )
npwqx = n_plane_waves( gcutw, q_grid%np, q_grid%p_cart, g, ngm )
!
ALLOCATE( gq2kin(npwqx) )
ALLOCATE( ngq(q_grid%nps) )
ALLOCATE( igq_q(npwqx,q_grid%nps) )
ALLOCATE( ngq(q_grid%np) )
ALLOCATE( igq_q(npwqx,q_grid%np) )
!ALLOCATE( igq_l2g(npwqx,q_grid%nps) )
igq_q(:,:) = 0
!igq_l2g(:,:) = 0
DO iq = 1, q_grid%nps
CALL gq_sort( q_grid%xp_cart(:,iq), ngm, g, gcutw, ngq(iq), igq_q(:,iq), gq2kin )
DO iq = 1, q_grid%np
CALL gq_sort( q_grid%p_cart(:,iq), ngm, g, gcutw, ngq(iq), igq_q(:,iq), gq2kin )
!CALL gq_l2gmap( ngm, ig_l2g(1), ngq(iq), igq_q(1,iq), igq_l2g(1,iq) )
ENDDO
!
......@@ -74,7 +74,7 @@ SUBROUTINE set_npwq()
!
! ... compute the global number of q+G vectors for each q-point
!
ALLOCATE( ngq_g(q_grid%nps) )
ALLOCATE( ngq_g(q_grid%np) )
!
ngq_g = 0
ngq_g(:) = ngq(:)
......@@ -86,7 +86,7 @@ SUBROUTINE set_npwq()
! ... compute the maximum G vector index among all q+G in processors
!
npwq_g = 0
DO iq = 1, q_grid%nps
DO iq = 1, q_grid%np
DO ig = 1, ngq(iq)
npwq_g = MAX( npwq_g, ig_l2g(igq_q(ig,iq)) )
ENDDO
......
......@@ -22,7 +22,90 @@ MODULE types_bz_grid
!
TYPE(bz_grid) :: k_grid
TYPE(bz_grid) :: q_grid
TYPE(bz_grid) :: kmq_grid
TYPE(bz_grid) :: kpq_grid
!TYPE(bz_grid) :: kmq_grid
!TYPE(bz_grid) :: kpq_grid
!
CONTAINS
!
!
FUNCTION findG(g0, unit_type) RETURN(ig0)
!
! ... ig0 is the index of G (unit_type = [ "cryst", "cart"])
! ... if on exit ig0 == 0 --> G is not found
!
USE cell_base, ONLY : bg
USE gvecs, ONLY : g, ngms
USE constants, ONLY : eps8
!
IMPLICIT NONE
!
! I/O
!
REAL(DP), INTENT(IN) :: g0(3)
CHARACTER(LEN=*), INTENT(IN) :: unit_type
!
! Workspace
!
REAL(DP) :: gtemp(3)
INTEGER :: ipol, ig, ig0
!
SELECT CASE(unit_type)
CASE("cryst","cart")
CASE DEFAULT
CALL errore( 1, "unit_type not supported, supported only cryst or cart" )
END SELECT
!
gtemp = g0
IF( unit_type == "cryst" ) CALL cryst_to_cart( 1, gtemp, bg, 1)
!
! gtemp is in cart
!
ig0 = 0
DO ig = 1, ngms
IF ( ALL ( ABS( g(:,ig) - g0(:) ) < eps8 ) ) THEN
ig0 = ig
ENDIF
ENDDO
!
END FUNCTION
!
!
SUBROUTINE compute_phase( g0, unit_type, phase )
!
! ... phase(r) = exp(-iG_0*r) (allocated externally)
!
USE gvecs, ONLY : nls
USE fft_base, ONLY : dffts
USE fft_interfaces, ONLY : invfft
!
IMPLICIT NONE
!
! I/O
!
REAL(DP), INTENT(IN) :: g0(3)
CHARACTER(LEN=*), INTENT(IN) :: unit_type
COMPLEX(DP), INTENT(OUT) :: phase(:)
!
! Workspace
!
INTEGER :: ipol, ig, ig0
!
SELECT CASE(unit_type)
CASE("cryst","cart")
CASE DEFAULT
CALL errore( 1, "unit_type not supported, supported only cryst or cart" )
END SELECT
!
ig0 = findG(g0,unit_type)
IF( ig0 == 0 ) CALL errore( 1, "G0 not found" )
!
phase = (0._DP, 0._DP)
phase( nls(ig0) ) = (1._DP, 0._DP)
!
! phase = exp(-iG_0*r)
CALL invfft( 'Wave', phase, dffts )
phase(1:dffts%nnr) = DCONJG( phase(1:dffts%nnr) )
!
END SUBROUTINE
!
END MODULE
......@@ -233,7 +233,7 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
USE io_push, ONLY : io_push_bar,io_push_value,io_push_title
USE distribution_center, ONLY : pert,ifr,rfr,aband
USE class_bz_grid, ONLY : bz_grid
USE types_bz_grid, ONLY : k_grid, q_grid, kmq_grid
USE types_bz_grid, ONLY : k_grid, q_grid
!
IMPLICIT NONE
!
......@@ -246,7 +246,7 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
!
! Workspace
!
INTEGER :: iks,iq,ikqs,ib,ifreq,glob_ifreq,il,im,glob_im,ip
INTEGER :: ik,iks,iq,ikqs,ib,ifreq,glob_ifreq,il,im,glob_im,ip
INTEGER :: nbndval
!
REAL(DP),EXTERNAL :: integrate_imfreq
......@@ -256,12 +256,11 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
!
COMPLEX(DP) :: partial_b,partial_h
REAL(DP) :: segno, enrg
REAL(DP) :: kmq(3), g0(3)
COMPLEX(DP) :: residues_b,residues_h
LOGICAL :: this_is_a_pole
LOGICAL :: l_gammaq
!
TYPE(bz_grid) :: k1_grid,q_grid_aux
!
! PRINT TITLE of CALC
!
IF(l_verbose) CALL io_push_title('Sigma_C')
......@@ -270,11 +269,6 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
!
CALL mp_barrier( world_comm )
!
! k1_grid = bz_grid()
! CALL k1_grid%init('K')
! q_grid_aux = bz_grid()
! CALL q_grid_aux%init_q( k_grid, k1_grid )
!
! ZERO
!
sigma_corr = 0._DP
......@@ -292,16 +286,19 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
! LOOP
!
DO iks = 1, k_grid%nps ! KPOINT-SPIN (MATRIX ELEMENT)
ik = k_grid%ip(iks)
!
DO ib = qp_bandrange(1), qp_bandrange(2)
!
partial_h = 0._DP
partial_b = 0._DP
!
DO iq = 1, q_grid%nps ! Q-POINT
DO iq = 1, q_grid%np ! Q-POINT
!
ikqs = kmq_grid%index_kq(iks,iq)
l_gammaq = q_grid%l_gammap(iq)
k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
ikqs = k_grid%find( kmq, 'cart' )
!ikqs = kmq_grid%index_kq(iks,iq)
l_gammaq = q_grid%l_pIsGamma(iq)
nbndval = nbnd_occ(ikqs)
!
! HEAD PART
......@@ -321,7 +318,7 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
DO im = 1, aband%nloc
glob_im = aband%l2g(im)
enrg = et(glob_im,ikqs) - energy(ib,iks)
partial_b = partial_b + z_body1_ifr_q(im,ifreq,ib,iks,iq)*integrate_imfreq(ifreq,enrg)*q_grid%wp(iq)
partial_b = partial_b + z_body1_ifr_q(im,ifreq,ib,iks,iq)*integrate_imfreq(ifreq,enrg)*q_grid%weight(iq)
ENDDO
ENDDO
!
......@@ -334,7 +331,7 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
DO ip = 1, pert%nloc
DO il = 1, n_lanczos
enrg = d_diago_q(il,ip,ib,iks,iq) - energy(ib,iks)
partial_b = partial_b + z_body2_ifr_q(il,ip,ifreq,ib,iks,iq)*integrate_imfreq(ifreq,enrg)*q_grid%wp(iq)
partial_b = partial_b + z_body2_ifr_q(il,ip,ifreq,ib,iks,iq)*integrate_imfreq(ifreq,enrg)*q_grid%weight(iq)
ENDDO
ENDDO
ENDDO
......@@ -371,6 +368,7 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
!
DO iks = 1, k_grid%nps
!
ik = k_grid%ip(iks)
!
DO ib = qp_bandrange(1), qp_bandrange(2)
!
......@@ -379,11 +377,13 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
residues_b = 0._DP
residues_h = 0._DP
!
DO iq = 1, q_grid%nps
DO iq = 1, q_grid%np
!
ikqs = kmq_grid%index_kq(iks,iq)
k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' )
ikqs = k_grid%find( kmq, 'cart' )
!ikqs = kmq_grid%index_kq(iks,iq)
l_gammaq = q_grid%l_pIsGamma(iq)
nbndval = nbnd_occ(ikqs)
l_gammaq = q_grid%l_gammap(iq)
!