TABLE OF CONTENTS


ABINIT/calc_coh_comp [ Functions ]

[ Top ] [ Functions ]

NAME

 calc_coh_comp

FUNCTION

  Calculates the COH-like contribution to the self-energy when
  the extrapolar technique and the closure relation is used to
  reduce the number of empty states to be summed over in the Green
  function entering the definition of the GW self-energy.

INPUTS

 iqibz=index of the irreducible q-point in the array qibz, point which is
  related by a symmetry operation to the point q summed over (see csigme).
  This index is also used to treat the integrable coulombian singularity at q=0
 ngfft(18)=contain all needed information about 3D FFT for GW wavefuntions,
  see ~abinit/doc/variables/vargs.htm#ngfft
 nsig_ab=Number of components in the self-energy operator (1 for collinear magnetism)
 npwc=number of plane waves in $\tilde epsilon^{-1}$
 nspinor=Number of spinorial components.
 i_sz=contribution arising from the integrable coulombian singularity at q==0
 (see csigme for the method used), note that in case of 3-D systems the factor
 4pi in the coulombian potential is included in the definition of i_sz
 gvec(3,npwc)=G vectors in reduced coordinates
 vc_sqrt(npwc)= square root of the coulombian matrix elements for this q-point
 botsq = Plasmon-pole parameters
 otq  = PPm parameters

OUTPUT

 sigcohme=partial contribution to the matrix element of $<jb k|\Sigma_{COH}| kb k>$
  coming from this single q-point for completeness trick

SIDE EFFECTS

PARENTS

      calc_sigc_me

CHILDREN

      wrtout

SOURCE

1390 subroutine calc_coh_comp(iqibz,i_sz,same_band,nspinor,nsig_ab,ediff,npwc,gvec,&
1391 &  ngfft,nfftot,wfg2_jk,vc_sqrt,botsq,otq,sigcohme)
1392 
1393 
1394 !This section has been created automatically by the script Abilint (TD).
1395 !Do not modify the following lines by hand.
1396 #undef ABI_FUNC
1397 #define ABI_FUNC 'calc_coh_comp'
1398 !End of the abilint section
1399 
1400  implicit none
1401 
1402 !Arguments ------------------------------------
1403 !scalars
1404  integer,intent(in) :: iqibz,npwc,nsig_ab,nspinor,nfftot
1405  real(dp),intent(in) :: i_sz,ediff
1406  logical,intent(in) :: same_band
1407 !arrays
1408  integer,intent(in) :: gvec(3,npwc),ngfft(18)
1409  complex(gwpc),intent(in) :: botsq(npwc,npwc),otq(npwc,npwc)
1410  complex(gwpc),intent(in) :: vc_sqrt(npwc)
1411  complex(gwpc),intent(in) :: wfg2_jk(nfftot*nsig_ab)
1412  complex(gwpc),intent(out) :: sigcohme(nsig_ab)
1413 
1414 !Local variables-------------------------------
1415 !scalars
1416  integer,save :: enough=0
1417  integer :: ig,ig4,ig4x,ig4y,ig4z,igp,igmin,ispinor,ngfft1,ngfft2,ngfft3,spad,outofbox
1418 !arrays
1419  integer :: g2mg1(3)
1420 
1421 ! *************************************************************************
1422 
1423  DBG_ENTER("COLL")
1424 
1425  ! === Treat the case q --> 0 adequately ===
1426  ! TODO Better treatment of wings
1427  igmin=1 ; if (iqibz==1) igmin=2
1428  !
1429  ! === Partial contribution to the matrix element of Sigma_c ===
1430  ! * For nspinor==2, the closure relation reads:
1431  !  $\sum_s \psi_a^*(1)\psi_b(2) = \delta_{ab} \delta(1-2)$
1432  !  where a,b are the spinor components. As a consequence, Sigma_{COH} is always
1433  !  diagonal in spin-space and only diagonal matrix elements have to be calculated.
1434  ! MG  TODO wfg2_jk should be calculated on an augmented FFT box to avoid spurious wrapping of G1-G2.
1435  !
1436  ngfft1 = ngfft(1); ngfft2 = ngfft(2); ngfft3 = ngfft(3)
1437  sigcohme(:) = czero_gw
1438 
1439  do ispinor=1,nspinor
1440   spad=(ispinor-1)*nfftot
1441   outofbox=0
1442 
1443    do igp=igmin,npwc
1444      do ig=igmin,npwc
1445 
1446       g2mg1 = gvec(:,igp)-gvec(:,ig)
1447       if (ANY(g2mg1(:)>ngfft(1:3)/2) .or. ANY(g2mg1(:)<-(ngfft(1:3)-1)/2)) then
1448         outofbox = outofbox+1; CYCLE
1449       end if
1450 
1451       ig4x=MODULO(g2mg1(1),ngfft1)
1452       ig4y=MODULO(g2mg1(2),ngfft2)
1453       ig4z=MODULO(g2mg1(3),ngfft3)
1454       ig4= 1+ig4x+ig4y*ngfft1+ig4z*ngfft1*ngfft2
1455 
1456       !MG where is neta here, ediff, otq might be close to zero depending on gwecomp
1457       sigcohme(ispinor) = sigcohme(ispinor) + &
1458 &       half*wfg2_jk(spad+ig4)*vc_sqrt(ig)*vc_sqrt(igp) * botsq(ig,igp) / ( otq(ig,igp) * ( ediff -otq(ig,igp) ) )
1459      end do
1460    end do
1461 
1462    if (iqibz==1.and.same_band) then
1463      sigcohme(ispinor) = sigcohme(ispinor) + half*wfg2_jk(spad+ig4)*i_sz*botsq(1,1) / ( otq(1,1) * (ediff -otq(1,1)) )
1464    end if
1465  end do !ispinor
1466 
1467  if (outofbox/=0) then
1468    enough=enough+1
1469    if (enough<=50) then
1470      MSG_WARNING(sjoin('Number of G1-G2 pairs outside the G-sphere for Wfns: ',itoa(outofbox)))
1471      if (enough==50) then
1472        call wrtout(std_out,' ========== Stop writing Warnings ==========')
1473      end if
1474    end if
1475  end if
1476 
1477  DBG_EXIT("COLL")
1478 
1479 end subroutine calc_coh_comp

ABINIT/calc_sig_ppm_comp [ Functions ]

[ Top ] [ Functions ]

NAME

 calc_sig_ppm_comp

FUNCTION

 Calculating contributions to self-energy operator using a plasmon-pole model

INPUTS

  nomega=number of frequencies to consider
  npwc= number of G vectors in the plasmon pole
  npwc1= 1 if ppmodel==3, =npwc if ppmodel== 4, 1 for all the other cases
  npwc2= 1 if ppmodel==3, =1    if ppmodel== 4, 1 for all the other cases
  npwx=number of G vectors in rhotwgp
  ppmodel=plasmon pole model
  theta_mu_minus_e0i= $\theta(\mu-\epsilon_{k-q,b1,s}), defines if the state is occupied or not
  zcut=small imaginary part to avoid the divergence. (see related input variable)
  omegame0i(nomega)=frequencies where evaluate \Sigma_c ($\omega$ - $\epsilon_i$
  otq(npwc,npwc2)=plasmon pole parameters for this q-point
  botsq(npwc,npwc1)=plasmon pole parameters for this q-point
  eig(npwc,npwc)=the eigvectors of the symmetrized inverse dielectric matrix for this q point
   (first index for G, second index for bands)
  rhotwgp(npwx)=oscillator matrix elements divided by |q+G| i.e
    $\frac{\langle b1 k-q s | e^{-i(q+G)r | b2 k s \rangle}{|q+G|}$

OUTPUT

  sigcme(nomega) (to be described), only relevant if ppm3 or ppm4

  ket(npwc,nomega):

  In case of ppmodel==1,2 it contains

   ket(G,omega) = Sum_G2       conjg(rhotw(G)) * Omega(G,G2) * rhotw(G2)
                          ---------------------------------------------------
                            2 omegatw(G,G2) (omega-E_i + omegatw(G,G2)(2f-1))

NOTES

 Taken from old routine

PARENTS

      calc_sigc_me

CHILDREN

SOURCE

1949 subroutine calc_sig_ppm_comp(npwc,nomega,rhotwgp,botsq,otq,omegame0i_io,zcut,theta_mu_minus_e0i,ket,ppmodel,npwx,npwc1,npwc2)
1950 
1951 
1952 !This section has been created automatically by the script Abilint (TD).
1953 !Do not modify the following lines by hand.
1954 #undef ABI_FUNC
1955 #define ABI_FUNC 'calc_sig_ppm_comp'
1956 !End of the abilint section
1957 
1958  implicit none
1959 
1960 !Arguments ------------------------------------
1961 !scalars
1962  integer,intent(in) :: nomega,npwc,npwc1,npwc2,npwx,ppmodel
1963  real(dp),intent(in) :: omegame0i_io,theta_mu_minus_e0i,zcut
1964 !arrays
1965  complex(gwpc),intent(in) :: botsq(npwc,npwc1),rhotwgp(npwx),otq(npwc,npwc2)
1966  complex(gwpc),intent(inout) :: ket(npwc,nomega)
1967 
1968 !Local variables-------------------------------
1969 !scalars
1970  integer :: ig,igp,io
1971  real(dp) :: den,otw,twofm1_zcut
1972  complex(gwpc) :: num,rhotwgdp_igp
1973  logical :: fully_occupied,totally_empty
1974  character(len=500) :: msg
1975 !arrays
1976  complex(gwpc),allocatable :: ket_comp(:)
1977 
1978 !*************************************************************************
1979 
1980  if (ppmodel/=1.and.ppmodel/=2) then
1981    write(msg,'(a,i0,a)')' The completeness trick cannot be used when ppmodel is ',ppmodel,' It should be set to 1 or 2. '
1982    MSG_ERROR(msg)
1983  end if
1984 
1985  ABI_ALLOCATE(ket_comp,(npwc))
1986  ket_comp(:)=0.d0
1987 
1988  fully_occupied=(abs(theta_mu_minus_e0i-1.)<0.001)
1989  totally_empty=(abs(theta_mu_minus_e0i)<0.001)
1990 
1991  if(.not.(totally_empty)) then ! not totally empty
1992    twofm1_zcut=zcut
1993    do igp=1,npwc
1994      rhotwgdp_igp=rhotwgp(igp)
1995      do ig=1,npwc
1996        otw=DBLE(otq(ig,igp)) ! in principle otw -> otw - ieta
1997        num = botsq(ig,igp)*rhotwgdp_igp
1998 
1999        den = omegame0i_io-otw
2000        if (den**2>zcut**2) then
2001          ket_comp(ig) = ket_comp(ig) - num/(den*otw)*theta_mu_minus_e0i
2002        end if
2003      end do !ig
2004    end do !igp
2005  end if ! not totally empty
2006 
2007  if(.not.(fully_occupied)) then ! not fully occupied
2008    twofm1_zcut=-zcut
2009 
2010    do igp=1,npwc
2011      rhotwgdp_igp=rhotwgp(igp)
2012      do ig=1,npwc
2013        otw=DBLE(otq(ig,igp)) ! in principle otw -> otw - ieta
2014        num = botsq(ig,igp)*rhotwgdp_igp
2015 
2016        den = omegame0i_io-otw
2017        if (den**2>zcut**2) then
2018          ket_comp(ig) = ket_comp(ig) - num/(den*otw)*(1.-theta_mu_minus_e0i)
2019        end if
2020      end do !ig
2021    end do !igp
2022  end if ! not fully occupied
2023 
2024  do io=1,nomega
2025    ket(:,io)=ket(:,io)+0.5*ket_comp(:)
2026  end do
2027 
2028  ABI_DEALLOCATE(ket_comp)
2029 
2030 end subroutine calc_sig_ppm_comp

ABINIT/calc_sigc_cd [ Functions ]

[ Top ] [ Functions ]

NAME

 calc_sigc_cd

FUNCTION

 Calculate contributions to the self-energy operator with the contour deformation method.

INPUTS

  nomega=Total number of frequencies where $\Sigma_c$ matrix elements are evaluated.
  nomegae=Number of frequencies where $\epsilon^{-1}$ has been evaluated.
  nomegaei=Number of imaginary frequencies for $\epsilon^{-1}$ (non zero).
  nomegaer=Number of real frequencies for $\epsilon^{-1}$
  npwc=Number of G vectors for the correlation part.
  npwx=Number of G vectors in rhotwgp for each spinorial component.
  nspinor=Number of spinorial components.
  theta_mu_minus_e0i=1 if e0i is occupied, 0 otherwise. Fractional occupancy in case of metals.
  omegame0i(nomega)=Contains $\omega-\epsilon_{k-q,b1,\sigma}$
  epsm1q(npwc,npwc,nomegae)=Symmetrized inverse dielectric matrix (exchange part is subtracted).
  omega(nomegae)=Set of frequencies for $\epsilon^{-1}$.
  rhotwgp(npwx*nspinor)=Matrix elements: $<k-q,b1,\sigma|e^{-i(q+G)r} |k,b2,\sigma>*vc_sqrt$

OUTPUT

 ket(npwc,nomega)=Contains \Sigma_c(\omega)|\phi> in reciprocal space.

SIDE EFFECTS

 npoles_missing=Incremented with the number of poles whose contribution has not been taken into account due to
  limited frequency mesh used for W.

PARENTS

      calc_sigc_me,m_screen

CHILDREN

      spline,splint,xgemm,xgemv

SOURCE

1518 subroutine calc_sigc_cd(npwc,npwx,nspinor,nomega,nomegae,nomegaer,nomegaei,rhotwgp,&
1519 & omega,epsm1q,omegame0i,theta_mu_minus_e0i,ket,plasmafreq,npoles_missing,calc_poles,method)
1520 
1521 
1522 !This section has been created automatically by the script Abilint (TD).
1523 !Do not modify the following lines by hand.
1524 #undef ABI_FUNC
1525 #define ABI_FUNC 'calc_sigc_cd'
1526 !End of the abilint section
1527 
1528  implicit none
1529 
1530 !Arguments ------------------------------------
1531 !scalars
1532  integer,intent(in) :: nomega,nomegae,nomegaei,nomegaer,npwc,npwx,nspinor
1533  integer,intent(inout) :: npoles_missing
1534  real(dp),intent(in) :: theta_mu_minus_e0i,plasmafreq
1535 !arrays
1536  real(dp),intent(in) :: omegame0i(nomega)
1537  complex(dpc),intent(in) :: omega(nomegae)
1538  complex(gwpc),intent(in) :: epsm1q(npwc,npwc,nomegae)
1539  complex(gwpc),intent(in) :: rhotwgp(npwx*nspinor)
1540  complex(gwpc),intent(inout) :: ket(nspinor*npwc,nomega)
1541  logical, intent(in), optional :: calc_poles(nomega)
1542  integer, intent(in), optional :: method
1543 
1544 !Local variables-------------------------------
1545 !scalars
1546  integer, parameter :: FABIEN=1,TRAPEZOID=2,NSPLINE=3
1547  integer :: ii,ig,io,ios,ispinor,spadc,spadx,my_err,ierr,GK_LEVEL,INTMETHOD
1548  integer :: i,j
1549  real(dp) :: rt_imag,rt_real,local_one,local_zero
1550  real(dp) :: intsign,temp1,temp2,temp3,temp4
1551  real(dp) :: alph,inv_alph,beta,alphsq,betasq,inv_beta
1552  real(dp) :: re_intG,re_intK,im_intG,im_intK,GKttab,tau,ttil
1553  real(dp) :: ref,imf,r,s,r2,s2
1554  complex(dpc) :: ct,domegaleft,domegaright
1555  complex(gwpc) :: fact
1556 !arrays
1557  real(dp) :: omegame0i_tmp(nomega),tmp_x(2),tmp_y(2)
1558  real(dp) :: left(nomega),right(nomega)
1559  real(dp) :: tbeta(nomega),tinv_beta(nomega),tbetasq(nomega)
1560  real(dp) :: atermr(nomega),aterml(nomega),logup(nomega),logdown(nomega)
1561  real(dp) :: rtmp_r(nomegaer),rtmp_i(nomegaer)
1562  real(dp) :: ftab(nomegaei+2),ftab2(nomegaei+2),xtab(nomegaei+2),y(3,nomegaei+2)
1563  real(dp) :: work(nomegaei+2),work2(nomegaei+2),y2(3,nomegaei+2)
1564  complex(dpc) :: omega_imag(nomegaei+1)
1565  complex(gwpc) :: epsrho(npwc,nomegae),epsrho_imag(npwc,nomegaei+1)
1566  complex(gwpc) :: tfone(npwc,nomegaei+1),tftwo(npwc,nomegaei+1)
1567  complex(gwpc) :: weight(nomegaei+1,nomega)
1568  complex(gwpc) :: weight2(nomegaei,nomega)
1569  logical :: my_calc_poles(nomega)
1570  real(dp), allocatable :: KronN(:),KronW(:),GaussW(:),fint(:),fint2(:)
1571 !*************************************************************************
1572 
1573  my_calc_poles=.TRUE.; my_err=0
1574 
1575  ! Set integration method for imaginary axis
1576  INTMETHOD = FABIEN
1577  if (present(method)) then
1578    if (method==1) INTMETHOD = FABIEN
1579    if (method==2) INTMETHOD = TRAPEZOID
1580    if (method>2) then
1581      INTMETHOD = NSPLINE
1582      if (method==3) then
1583        GK_LEVEL = 15
1584        ABI_ALLOCATE(KronN,(GK_LEVEL))
1585        ABI_ALLOCATE(KronW,(GK_LEVEL))
1586        ABI_ALLOCATE(GaussW,(GK_LEVEL-8))
1587        ABI_ALLOCATE(fint,(GK_LEVEL))
1588        ABI_ALLOCATE(fint2,(GK_LEVEL))
1589        KronN(:) = Kron15N(:); KronW(:) = Kron15W(:); GaussW(:) = Gau7W(:)
1590      else if (method==4) then
1591        GK_LEVEL = 23
1592        ABI_ALLOCATE(KronN,(GK_LEVEL))
1593        ABI_ALLOCATE(KronW,(GK_LEVEL))
1594        ABI_ALLOCATE(GaussW,(GK_LEVEL-12))
1595        ABI_ALLOCATE(fint,(GK_LEVEL))
1596        ABI_ALLOCATE(fint2,(GK_LEVEL))
1597        KronN(:) = Kron23N(:); KronW(:) = Kron23W(:); GaussW(:) = Gau11W(:)
1598      else if (method>4) then
1599        GK_LEVEL = 31
1600        ABI_ALLOCATE(KronN,(GK_LEVEL))
1601        ABI_ALLOCATE(KronW,(GK_LEVEL))
1602        ABI_ALLOCATE(GaussW,(GK_LEVEL-16))
1603        ABI_ALLOCATE(fint,(GK_LEVEL))
1604        ABI_ALLOCATE(fint2,(GK_LEVEL))
1605        KronN(:) = Kron31N(:); KronW(:) = Kron31W(:); GaussW(:) = Gau15W(:)
1606      end if
1607    end if
1608  end if
1609 
1610  ! Avoid divergences in $\omega - \omega_s$.
1611  omegame0i_tmp(:)=omegame0i(:)
1612  do ios=1,nomega
1613    if (ABS(omegame0i_tmp(ios))<tol6) omegame0i_tmp(ios)=sign(tol6,omegame0i_tmp(ios))
1614  end do
1615 
1616  do ispinor=1,nspinor
1617    spadx=(ispinor-1)*npwx; spadc=(ispinor-1)*npwc
1618 
1619    ! Calculate $ \sum_{Gp} (\epsilon^{-1}_{G Gp}(\omega)-\delta_{G Gp}) \rhotwgp(Gp) $
1620 !$omp parallel do
1621    do io=1,nomegae
1622      call XGEMV('N',npwc,npwc,cone_gw,epsm1q(:,:,io),npwc,rhotwgp(1+spadx:),1,czero_gw,epsrho(:,io),1)
1623    end do
1624 
1625    ! Integrand along the imaginary axis.
1626    epsrho_imag(:,1)=epsrho(:,1)
1627    epsrho_imag(:,2:nomegaei+1)=epsrho(:,nomegaer+1:nomegae)
1628 
1629    ! Frequency mesh for integral along the imaginary axis.
1630    omega_imag(1)=omega(1)
1631    omega_imag(2:nomegaei+1)=omega(nomegaer+1:nomegae)
1632 
1633    ! Original implementation -- saved here for reference during development
1634    ! === Perform integration along the imaginary axis ===
1635    !do io=1,nomegaei+1
1636    !  if (io==1) then
1637    !    domegaleft  = omega_imag(io)
1638    !    domegaright =(omega_imag(io+1)-omega_imag(io  ))*half
1639    !  else if (io==nomegaei+1) then
1640    !    domegaleft  =(omega_imag(io  )-omega_imag(io-1))*half
1641    !    domegaright =(omega_imag(io  )-omega_imag(io-1))*half
1642    !  else
1643    !    domegaleft  =(omega_imag(io  )-omega_imag(io-1))*half
1644    !    domegaright =(omega_imag(io+1)-omega_imag(io  ))*half
1645    !  end if
1646    !  do ios=1,nomega
1647    !    omg2 = -AIMAG(omega_imag(io)+domegaright)/REAL(omegame0i_tmp(ios))
1648    !    omg1 = -AIMAG(omega_imag(io)-domegaleft )/REAL(omegame0i_tmp(ios))
1649    !    fact = ATAN(omg2)-ATAN(omg1)
1650    !    ket(spadc+1:spadc+npwc,ios)=ket(spadc+1:spadc+npwc,ios)+epsrho_imag(:,io)*fact
1651    !  end do
1652    !end do !io
1653 
1654    !ket(spadc+1:spadc+npwc,:)=ket(spadc+1:spadc+npwc,:)/pi
1655    ! ---------------- end of original implementation -----------------------
1656 
1657    select case(INTMETHOD)
1658    case(FABIEN)
1659      ! Hopefully more effective implementation MS 04.08.2011
1660      ! Perform integration along imaginary axis using BLAS
1661      ! First calculate first and last weights
1662      weight(1,:) = ATAN(-half*AIMAG(omega_imag(2))/REAL(omegame0i_tmp(:)))
1663      domegaleft  = (three*omega_imag(nomegaei+1)-omega_imag(nomegaei))
1664      domegaright = (omega_imag(nomegaei+1)+omega_imag(nomegaei))
1665      right(:)    = -AIMAG(omega_imag(nomegaei+1)-omega_imag(nomegaei))*REAL(omegame0i_tmp(:))
1666      left(:)     = quarter*AIMAG(domegaleft)*AIMAG(domegaright) &
1667 &                   +REAL(omegame0i_tmp(:))*REAL(omegame0i_tmp(:))
1668      weight(nomegaei+1,:) = ATAN(right(:)/left(:))
1669      ! Calculate the rest of the weights
1670      do io=2,nomegaei
1671        domegaleft  = (omega_imag(io  )+omega_imag(io-1))
1672        domegaright = (omega_imag(io+1)+omega_imag(io  ))
1673        right(:)    = -half*AIMAG(omega_imag(io+1)-omega_imag(io-1))*REAL(omegame0i_tmp(:))
1674        left(:)     = REAL(omegame0i_tmp(:))*REAL(omegame0i_tmp(:)) &
1675 &       +quarter*AIMAG(domegaleft)*AIMAG(domegaright)
1676        weight(io,:) = ATAN(right(:)/left(:))
1677      end do
1678 
1679      ! Use BLAS call to perform matrix-matrix multiplication and accumulation
1680      fact = CMPLX(piinv,zero)
1681 
1682      call xgemm('N','N',npwc,nomega,nomegaei+1,fact,epsrho_imag,npwc,&
1683 &     weight,nomegaei+1,cone_gw,ket(spadc+1:spadc+npwc,:),npwc)
1684 
1685    case(TRAPEZOID)
1686 !   Trapezoidal rule Transform omega coordinates
1687      alph     = plasmafreq
1688      alphsq   = alph*alph
1689      inv_alph = one/alph
1690 
1691      xtab(1:nomegaei+1) = AIMAG(omega_imag(:))/(AIMAG(omega_imag(:)) + alph)
1692      xtab(nomegaei+2)   = one
1693 
1694 !   Efficient trapezoidal rule with BLAS calls
1695      tbeta(:)     = REAL(omegame0i_tmp(:))
1696      tbetasq(:)   = tbeta(:)*tbeta(:)
1697      tinv_beta(:) = one/tbeta(:)
1698 
1699      do io=1,nomegaei
1700        atermr(:)    = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(io+1)-tbetasq(:))
1701        aterml(:)    = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(io  )-tbetasq(:))
1702        right(:)     = ATAN((atermr(:)-aterml(:))/(one+atermr(:)*aterml(:)))
1703        logup(:)     = ABS(((alphsq+tbetasq(:))*xtab(io+1)-two*tbetasq(:)) &
1704 &                     *xtab(io+1)+tbetasq(:))
1705        logdown(:)   = ABS(((alphsq+tbetasq(:))*xtab(io  )-two*tbetasq(:)) &
1706 &                     *xtab(io  )+tbetasq(:))
1707        ! Trapezoid integration weights
1708        weight(io,:)  = CMPLX(-(half*alph*tbeta(:)*LOG(logup(:)/logdown(:)) + tbetasq(:) &
1709 &                         *right(:))/(alphsq+tbetasq(:)),zero)
1710        weight2(io,:) = CMPLX(-right(:),zero)
1711        ! Linear interpolation coefficients for each section (sum over ig)
1712        tfone(:,io)   = (epsrho_imag(:,io+1)-epsrho_imag(:,io)) &
1713 &                      /(xtab(io+1)-xtab(io))
1714        tftwo(:,io)   = epsrho_imag(:,io) - tfone(:,io)*xtab(io)
1715      end do
1716 
1717      ! Calculate weights for asymptotic behaviour
1718      atermr(:)   = alph*tinv_beta(:)
1719      aterml(:)   = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(nomegaei+1)-tbetasq(:))
1720      logup(:)    = alphsq*xtab(nomegaei+1)*xtab(nomegaei+1)
1721      logdown(:)  = ABS(((alphsq+tbetasq(:))*xtab(nomegaei+1)-two*tbetasq(:)) &
1722 &                   *xtab(nomegaei+1)+tbetasq(:))
1723      right(:)     = ATAN((atermr(:)-aterml(:))/(one+atermr(:)*aterml(:)))
1724      weight (nomegaei+1,:) = CMPLX(-(half*(alphsq*tinv_beta(:)*LOG(logdown(:)/logup(:)) &
1725 &     - tbeta(:)*LOG(xtab(nomegaei+1)*xtab(nomegaei+1))) - alph*right(:)),zero)
1726      tfone(:,nomegaei+1) = -(zero-epsrho_imag(:,nomegaei+1)*AIMAG(omega_imag(nomegaei+1))) &
1727                            /(one-xtab(nomegaei+1))
1728 
1729      ! Use BLAS call to perform matrix-matrix multiplication and accumulation
1730      fact = CMPLX(piinv,zero)
1731 
1732      call xgemm('N','N',npwc,nomega,nomegaei+1,fact,tfone,npwc,&
1733 &     weight ,nomegaei+1,cone_gw,ket(spadc+1:spadc+npwc,:),npwc)
1734      call xgemm('N','N',npwc,nomega,nomegaei  ,fact,tftwo,npwc,&
1735 &     weight2,nomegaei  ,cone_gw,ket(spadc+1:spadc+npwc,:),npwc)
1736 
1737    case(NSPLINE)
1738      ! Natural spline followed by Gauss-Kronrod
1739      ! Transform omega coordinates
1740      alph     = plasmafreq
1741      alphsq   = alph*alph
1742      inv_alph = one/alph
1743 
1744      xtab(1:nomegaei+1) = AIMAG(omega_imag(:))/(AIMAG(omega_imag(:)) + alph)
1745      xtab(nomegaei+2)   = one
1746 
1747 ! Gauss-Kronrod integration of spline fit of f(t)/(1-t) in transformed space
1748 ! *** OPENMP SECTION *** Added by MS
1749 !!$OMP PARALLEL DO PRIVATE(ig,ftab,ftab2,s,s2,r,r2,y,y2,work,work2,beta,betasq,inv_beta, &
1750 !!$OMP  intsign,io,ii,i,j,re_intG,re_intK,im_intG,im_intK,temp1,temp2,temp3,temp4, &
1751 !!$OMP  ttil,tau,ref,fint,imf,fint2,GKttab)
1752      do ig=1,npwc
1753        ! Spline fit
1754        ftab (1:nomegaei+1) =  REAL(epsrho_imag(ig,1:nomegaei+1))/(one-xtab(1:nomegaei+1))
1755        ftab2(1:nomegaei+1) = AIMAG(epsrho_imag(ig,1:nomegaei+1))/(one-xtab(1:nomegaei+1))
1756        ftab (nomegaei+2)   = zero; ftab2(nomegaei+2) = zero
1757        ! Explicit calculation of spline coefficients
1758        s  = zero; s2 = zero
1759        do i = 1, nomegaei+2-1
1760          r  = ( ftab (i+1) - ftab (i) ) / ( xtab(i+1) - xtab(i) )
1761          r2 = ( ftab2(i+1) - ftab2(i) ) / ( xtab(i+1) - xtab(i) )
1762          y (2,i) = r  - s; y2(2,i) = r2 - s2
1763          s  = r; s2 = r2
1764        end do
1765        s = zero; s2 = zero
1766        r = zero; r2 = zero
1767        y(2,1) = zero; y2(2,1) = zero
1768        y(2,nomegaei+2) = zero; y2(2,nomegaei+2) = zero
1769        do i = 2, nomegaei+2-1
1770          y (2,i) = y (2,i) + r  * y (2,i-1)
1771          y2(2,i) = y2(2,i) + r2 * y2(2,i-1)
1772          work (i) = two * ( xtab(i-1) - xtab(i+1) ) - r  * s
1773          work2(i) = two * ( xtab(i-1) - xtab(i+1) ) - r2 * s2
1774          s = xtab(i+1) - xtab(i)
1775          s2 = s
1776          r  = s  / work (i)
1777          r2 = s2 / work2(i)
1778        end do
1779        do j = 2, nomegaei+2-1
1780          i = nomegaei+2+1-j
1781          y (2,i) = ( ( xtab(i+1) - xtab(i) ) * y (2,i+1) - y (2,i) ) / work (i)
1782          y2(2,i) = ( ( xtab(i+1) - xtab(i) ) * y2(2,i+1) - y2(2,i) ) / work2(i)
1783        end do
1784        do i = 1, nomegaei+2-1
1785          s = xtab(i+1) - xtab(i); s2 = s;
1786          r = y(2,i+1) - y(2,i); r2 = y2(2,i+1) - y2(2,i);
1787          y(3,i) = r / s; y2(3,i) = r2 / s2;
1788          y(2,i) = three * y(2,i); y2(2,i) = three * y2(2,i);
1789          y (1,i) = ( ftab (i+1) - ftab (i) ) / s  - ( y (2,i) + r  ) * s
1790          y2(1,i) = ( ftab2(i+1) - ftab2(i) ) / s2 - ( y2(2,i) + r2 ) * s2
1791        end do
1792        ! End of spline interpolation
1793        do ios=1,nomega
1794          beta     = REAL(omegame0i_tmp(ios))
1795          betasq   = beta*beta
1796          inv_beta = one/beta
1797          intsign = sign(half*piinv,beta)
1798          beta = ABS(beta)
1799          io = 1; re_intG = zero; re_intK = zero; im_intG = zero; im_intK = zero
1800          do ii=1,GK_LEVEL
1801            do
1802              GKttab = two*alph*xtab(io+1)/(beta-(beta-alph)*xtab(io+1))-one
1803              if (GKttab > KronN(ii)) EXIT
1804              io = io + 1
1805            end do
1806            temp1     = half*(KronN(ii)+one)
1807            temp2     = temp1 - half
1808            temp3     = temp2*temp2
1809            temp4     = half/(temp3 + quarter)
1810            ttil      = beta*temp1/(alph-(alph-beta)*temp1)
1811            tau       = ttil - xtab(io)
1812            ref       = ftab (io) + tau*(y (1,io)+tau*(y (2,io)+tau*y (3,io)))
1813            fint (ii) = -ref*(one-ttil)*temp4
1814            imf       = ftab2(io) + tau*(y2(1,io)+tau*(y2(2,io)+tau*y2(3,io)))
1815            fint2(ii) = -imf*(one-ttil)*temp4
1816            re_intK   = KronW(ii)*fint (ii)
1817            im_intK   = KronW(ii)*fint2(ii)
1818            ket(spadc+ig,ios) = ket(spadc+ig,ios)+intsign*CMPLX(re_intK,im_intK)
1819            end do ! ii
1820        end do !ios
1821      end do !ig
1822 !!$OMP END PARALLEL DO
1823 
1824    end select !intmethod
1825 
1826    local_one = one
1827    local_zero = zero
1828 
1829    ! ============================================
1830    ! ==== Add contribution coming from poles ====
1831    ! ============================================
1832    ! First see if the contribution has been checked before the routine is entered
1833    if (present(calc_poles)) then
1834      my_calc_poles = calc_poles
1835    else ! Otherwise check locally if there is a contribution
1836      do ios=1,nomega
1837        if (omegame0i_tmp(ios)>tol12) then
1838          if ((local_one-theta_mu_minus_e0i)<tol12) my_calc_poles(ios) = .FALSE.
1839        end if
1840        if (omegame0i_tmp(ios)<-tol12) then
1841          if (theta_mu_minus_e0i<tol12) my_calc_poles(ios) = .FALSE.
1842        end if
1843      end do !ios
1844    end if
1845 
1846    if (ANY(my_calc_poles(:))) then ! Make sure we only enter if necessary
1847 ! *** OPENMP SECTION *** Added by MS
1848 !!OMP !write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',xomp_get_num_threads()
1849 !$OMP PARALLEL SHARED(npwc,nomega,nomegaer,theta_mu_minus_e0i,spadc,local_one,local_zero, &
1850 !$OMP                    omega,epsrho,omegame0i_tmp,ket,my_calc_poles) &
1851 !$OMP PRIVATE(ig,ios,rtmp_r,rtmp_i,tmp_x,tmp_y,rt_real,rt_imag,ct,ierr) REDUCTION(+:my_err)
1852 !!OMP $ write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',xomp_get_num_threads()
1853 !$OMP DO
1854      do ig=1,npwc
1855        ! Prepare the spline interpolation by filling at once the arrays rtmp_r, rtmp_i
1856        call spline(DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),nomegaer,local_zero,local_zero,rtmp_r)
1857        call spline(DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),nomegaer,local_zero,local_zero,rtmp_i)
1858        ! call spline_complex( DBLE(omega(1:nomegaer)), epsrho(ig,1:nomegaer), nomegaer, zero, zero, rtmp )
1859 
1860        do ios=1,nomega
1861          if (.NOT.my_calc_poles(ios)) CYCLE
1862 
1863          ! Interpolate real and imaginary part of epsrho at |omegame0i_tmp|.
1864          tmp_x(1) = ABS(omegame0i_tmp(ios))
1865          call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),rtmp_r,1,tmp_x,tmp_y,ierr=ierr)
1866          if (ig==1.and.ispinor==1) my_err = my_err + ierr
1867          rt_real = tmp_y(1)
1868 
1869          tmp_x(1) = ABS(omegame0i_tmp(ios))
1870          call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),rtmp_i,1,tmp_x,tmp_y)
1871          rt_imag = tmp_y(1)
1872          !!call splint_complex(nomegaer,DBLE(omega(1:nomegaer)),epsrho(ig,1:nomegaer),rtmp,1,tmp_x,yfit)
1873 
1874          ct=DCMPLX(rt_real,rt_imag)
1875 
1876          if (omegame0i_tmp(ios)>tol12) then
1877            ket(spadc+ig,ios)=ket(spadc+ig,ios)+ct*(local_one-theta_mu_minus_e0i)
1878          end if
1879          if (omegame0i_tmp(ios)<-tol12) then
1880            ket(spadc+ig,ios)=ket(spadc+ig,ios)-ct*theta_mu_minus_e0i
1881          end if
1882 
1883        end do !ios
1884      end do !ig
1885 !$OMP END DO
1886 !$OMP END PARALLEL
1887    end if ! ANY(my_calc_poles)
1888  end do !ispinor
1889 
1890  npoles_missing = npoles_missing + my_err
1891 
1892  if (INTMETHOD>2) then
1893    ABI_DEALLOCATE(KronN)
1894    ABI_DEALLOCATE(KronW)
1895    ABI_DEALLOCATE(GaussW)
1896    ABI_DEALLOCATE(fint)
1897    ABI_DEALLOCATE(fint2)
1898  end if
1899 
1900 end subroutine calc_sigc_cd

ABINIT/calc_sigc_me [ Functions ]

[ Top ] [ Functions ]

NAME

 calc_sigc_me

FUNCTION

 Calculate diagonal and off-diagonal matrix elements of the self-energy operator.

INPUTS

 sigmak_ibz=Index of the k-point in the IBZ.
 minbnd, maxbnd= min and Max band index for GW correction (for this k-point)
 Dtset <type(dataset_type)>=all input variables in this dataset
    %iomode
    %paral_kgb
    %nspinor=Number of spinorial components.
    %gwcomp=If 1 use an extrapolar approximation to accelerate convergence.
    %gwencomp=Extrapolar energy.
 Er <Epsilonm1_results> (see the definition of this structured datatype)
    %mqmem=if 0 use out-of-core method in which a single q-slice of espilon is read inside the loop over k
      much slower but it requires less memory
    %nomega_i=Number of imaginary frequencies.
    %nomega_r=Number of real frequencies.
    %nomega=Total number of frequencies.
 Gsph_c<gsphere_t>= info on G-sphere for Sigma_c
 Gsph_Max<gsphere_t>= info on biggest G-sphere
    %nsym=number of symmetry operations
    %rottb(ng,timrev,nsym)=index of (IS) G where I is the identity or the inversion
      operation and G is one of the ng vectors in reciprocal space
    %timrev=2 if time-reversal symmetry is used, 1 otherwise
    %gvec(3,ng)=integer coordinates of each plane wave in reciprocal space
 ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated
 Ltg_k datatype containing information on the little group
 Kmesh <kmesh_t>
    %nbz=Number of points in the BZ
    %nibz=Number of points in IBZ
    %kibz(3,nibz)=k-point coordinates, irreducible Brillouin zone
    %kbz(3,nbz)=k-point coordinates, full Brillouin zone
    %ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
    %ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered
    %ktabp(nbz)= phase factor associated to tnons
 gwc_ngfft(18)=Information about 3D FFT for the oscillator strengths used for the correlation part,
 Vcp <vcoul_t datatype> containing information on the cutoff technique
    %vc_sqrt(npwc,nqibz)= square-root of the coulombian potential for q-points in the IBZ
 Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
 Pawang <type(pawang_type)>=paw angular mesh and related data
 Psps <type(pseudopotential_type)>=variables related to pseudopotentials
    %usepaw=1 for PAW, 0 for NC pseudopotentials.
 Qmesh <kmesh_t> : datatype gathering information of the q-mesh used
    %ibz=q points where $\tilde\epsilon^{-1}$ has been computed
    %bz(3,nqbz)=coordinates of all q-points in BZ
 Sigp <sigparams_t> (see the definition of this structured datatype)
 Cryst<crystal_t>=Info on unit cell and symmetries
    %natom=number of atoms in unit cell
    %ucvol=unit cell volume
    %nsym=number of symmetry operations
    %typat(natom)=type of each atom
 PPm<ppmodel_t>= Datatype gathering information on the Plasmonpole technique (see also ppm_get_qbz).
    %model=type of ppmodel
    %npwc=number of G-vectors for the correlation part.
    %dm2_otq =size of second dimension of otq array
    %dm2_bots=size of second dimension of botsq arrays
    %dm_eig  =size of second dimension of eig arrays
 QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
  eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
  occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
  Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
 allQP_sym(Wfd%nkibz,Wfd%nsppol)<esymm_t>=Datatype collecting data on the irreducible representaions of the
    little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states.
  Sr=sigma_t (see the definition of this structured datatype)
  use_aerhor=1 is aepaw_rhor is used, 0 otherwise.
  aepaw_rhor(rho_nfftot,Wfd%nspden*use_aerhor)=AE PAW density used to generate PPmodel paramenters if mqmem==0

OUTPUT

NOTES

  1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) [[cite:Gigy1986]] is included.
  2) The calculation of energy derivative is based on finite elements.
  3) On the symmetrization of Sigma matrix elements ***/
        If  Sk = k+G0 then  M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G}   (k,q)
        If -Sk = k+G0 then  M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q)

     Notice the absence of G0 in the expression. Moreover, when we sum over the little group, it turns out
     that there is a cancellation of the phase factor associated to the non-symmorphic operations due to a
     similar term coming from the symmetrization of \epsilon^{-1}. Mind however that the nonsymmorphic phase
     has to be considered when epsilon^-1 is reconstructed starting from the q-points in the IBZ.

  4) The unitary transformation relating wavefunctions
     at symmetric k-points should be taken into account during the symmetrization
     of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however,
     it is possible to make an invariant by just including all the degenerate states and
     averaging the final results over the degenerate subset.

PARENTS

      sigma

CHILDREN

      calc_coh_comp,calc_sig_ppm,calc_sig_ppm_comp,calc_sigc_cd,calc_wfwfg
      coeffs_gausslegint,cwtime,epsm1_symmetrizer,epsm1_symmetrizer_inplace
      esymm_symmetrize_mels,findqg0,get_bz_item,get_epsm1,get_gftt
      gsph_fft_tabs,littlegroup_print,paw_cross_rho_tw_g,paw_rho_tw_g
      paw_symcprj,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawpwij_free
      pawpwij_init,ppm_get_qbz,rho_tw_g,rotate_fft_mesh,setup_ppmodel
      sigma_distribute_bks,timab,wfd_change_ngfft,wfd_get_cprj
      wfd_get_many_ur,wfd_get_ur,wfd_paw_get_aeur,wrtout,xmpi_max,xmpi_sum

SOURCE

 178 subroutine calc_sigc_me(sigmak_ibz,ikcalc,nomega_sigc,minbnd,maxbnd,&
 179 & Dtset,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_Max,Gsph_c,Vcp,Kmesh,Qmesh,Ltg_k,&
 180 & PPm,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,allQP_sym,&
 181 & gwc_ngfft,rho_ngfft,rho_nfftot,rhor,use_aerhor,aepaw_rhor,sigcme_tmp)
 182 
 183 
 184 !This section has been created automatically by the script Abilint (TD).
 185 !Do not modify the following lines by hand.
 186 #undef ABI_FUNC
 187 #define ABI_FUNC 'calc_sigc_me'
 188 !End of the abilint section
 189 
 190  implicit none
 191 
 192 !Arguments ------------------------------------
 193 !scalars
 194  integer,intent(in) :: sigmak_ibz,ikcalc,rho_nfftot,nomega_sigc,minbnd,maxbnd
 195  integer,intent(in) :: use_aerhor
 196  type(crystal_t),intent(in) :: Cryst
 197  type(ebands_t),target,intent(in) :: QP_BSt
 198  type(kmesh_t),intent(in) :: Kmesh,Qmesh
 199  type(vcoul_t),intent(in) :: Vcp
 200  type(dataset_type),intent(in) :: Dtset
 201  type(epsilonm1_results),intent(inout) :: Er
 202  type(gsphere_t),intent(in) :: Gsph_Max,Gsph_c
 203  type(littlegroup_t),intent(in) :: Ltg_k
 204  type(ppmodel_t),intent(inout) :: PPm
 205  type(Pseudopotential_type),intent(in) :: Psps
 206  type(pawang_type),intent(in) :: pawang
 207  type(sigparams_t),target,intent(in) :: Sigp
 208  type(sigma_t),intent(in) :: Sr
 209  type(wfd_t),target,intent(inout) :: Wfd,Wfdf
 210 !arrays
 211  integer,intent(in) :: gwc_ngfft(18),rho_ngfft(18)
 212  real(dp),intent(in) :: rhor(rho_nfftot,Wfd%nspden)
 213  real(dp),intent(in) :: aepaw_rhor(rho_nfftot,Wfd%nspden*use_aerhor)
 214  complex(dpc),intent(out) :: sigcme_tmp(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Wfd%nsppol*Sigp%nsig_ab)
 215  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
 216  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
 217  type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol)
 218  type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw)
 219  type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom*Psps%usepaw)
 220 
 221 !Local variables ------------------------------
 222 !scalars
 223  integer,parameter :: tim_fourdp2=2,ndat1=1
 224  integer :: npw_k,iab,ib,ib1,ib2,ierr,ig,iggp,igp,ii,iik,itim_q,i1,i2,npls
 225  integer :: ik_bz,ik_ibz,io,iiw,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx
 226  integer :: band,band1,band2,idle,rank,jik,jk_bz,jk_ibz,kb,nspinor
 227  integer :: nomega_tot,nq_summed,ispinor,ibsp,dimcprj_gw,npwc
 228  integer :: spad,spadc,spadc1,spadc2,irow,my_nbks,ndegs,wtqm,wtqp,mod10
 229  integer :: isym_kgw,isym_ki,gwc_mgfft,use_padfft,gwc_fftalga,gwc_nfftot,nfftf,mgfftf,use_padfftf
 230  integer :: iwc,ifft
 231  real(dp) :: cpu_time,wall_time,gflops
 232  real(dp) :: e0i,fact_sp,theta_mu_minus_e0i,tol_empty,z2,en_high,gw_gsq,w_localmax,w_max
 233  complex(dpc) :: ctmp,omegame0i2_ac,omegame0i_ac,ph_mkgwt,ph_mkt
 234  logical :: iscompatibleFFT,q_is_gamma
 235  character(len=500) :: msg,sigma_type
 236  complex(gwpc),allocatable :: botsq(:,:),otq(:,:),eig(:,:)
 237 !arrays
 238  integer :: g0(3),spinor_padc(2,4),got(Wfd%nproc)
 239  integer,allocatable :: proc_distrb(:,:,:),extrapolar_distrb(:,:,:,:),degtab(:,:,:)
 240  integer,allocatable :: igfftcg0(:),gw_gfft(:,:),gw_gbound(:,:),irottb(:,:),ktabr(:,:)
 241  integer,allocatable :: igfftfcg0(:),gboundf(:,:),ktabrf(:,:),npoles_missing(:)
 242  real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),omegap(Er%nomega_i),omegap2(Er%nomega_i),q0(3),tsec(2),qbz(3)
 243  real(dp) :: gl_knots(Er%nomega_i),gl_wts(Er%nomega_i)
 244  real(dp) :: spinrot_kbz(4),spinrot_kgw(4)
 245  real(dp),ABI_CONTIGUOUS pointer :: qp_ene(:,:,:),qp_occ(:,:,:)
 246  real(dp),allocatable :: omegame0i(:), w_maxval(:)
 247  complex(gwpc) :: sigcohme(Sigp%nsig_ab)
 248  complex(gwpc),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:)
 249  complex(gwpc),allocatable :: botsq_conjg_transp(:,:),ac_epsm1cqwz2(:,:,:)
 250  complex(gwpc),allocatable :: epsm1_qbz(:,:,:),epsm1_trcc_qbz(:,:,:), epsm1_tmp(:,:)
 251  complex(gwpc),allocatable :: ac_integr(:,:,:),sigc_ket(:,:),ket1(:,:),ket2(:,:)
 252  complex(gwpc),allocatable :: herm_sigc_ket(:,:),aherm_sigc_ket(:,:)
 253  complex(gwpc),allocatable :: rhotwg_ki(:,:)
 254  complex(gwpc),allocatable :: sigcme2(:,:),sigcme_3(:),sigcme_new(:),sigctmp(:,:)
 255  complex(gwpc),allocatable :: wfr_bdgw(:,:),ur_ibz(:),wf1swf2_g(:),usr_bz(:)
 256  complex(gwpc),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:)
 257  complex(gwpc),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:)
 258  complex(gwpc),allocatable :: otq_transp(:,:)
 259  complex(gwpc),ABI_CONTIGUOUS pointer :: cg_jb(:),cg_sum(:)
 260  complex(dpc),allocatable :: sym_cme(:,:,:,:),sigc(:,:,:,:,:)
 261  logical :: rank_mask(Wfd%nproc),can_symmetrize(Wfd%nsppol)
 262 !logical :: me_calc_poles(Sr%nomega_r+Sr%nomega4sd)
 263  type(sigijtab_t),pointer :: Sigcij_tab(:)
 264  type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:)
 265  type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:)
 266  type(esymm_t),pointer :: QP_sym(:)
 267 
 268 !************************************************************************
 269 
 270  DBG_ENTER("COLL")
 271 
 272  ! Initial check
 273  ABI_CHECK(Sr%nomega_r==Sigp%nomegasr,"")
 274  ABI_CHECK(Sr%nomega4sd==Sigp%nomegasrd,"")
 275  ABI_CHECK(Sigp%npwc==Gsph_c%ng,"")
 276  ABI_CHECK(Sigp%npwvec==Gsph_Max%ng,"")
 277 
 278  call timab(424,1,tsec) ! calc_sigc_me
 279  call timab(431,1,tsec) ! calc_sigc_me
 280  call timab(432,1,tsec) ! Init
 281  call cwtime(cpu_time,wall_time,gflops,"start")
 282 
 283  ! Initialize MPI variables
 284  qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ
 285 
 286  ! Extract the symmetries of the bands for this k-point
 287  QP_sym => allQP_sym(sigmak_ibz,1:Wfd%nsppol)
 288 
 289  ! Index of the GW point in the BZ array, its image in IBZ and time-reversal ===
 290  jk_bz=Sigp%kptgw2bz(ikcalc)
 291  call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
 292  !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
 293 
 294  ! TODO: the new version based of get_uug won't suppporte kptgw vectors that are not in
 295  ! the IBZ since one should perform the rotation before entering the band loop
 296  ! In the old version, the rotation was done in rho_tw_g
 297  !ABI_CHECK(jik==1,"jik!=1")
 298  !ABI_CHECK(isym_kgw==1,"isym_kgw!=1")
 299  !ABI_CHECK((ABS(ph_mkgwt - cone) < tol12),"ph_mkgwt!")
 300 
 301  spinrot_kgw=Cryst%spinrot(:,isym_kgw)
 302  ib1=minbnd; ib2=maxbnd
 303 
 304  write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,&
 305 &  ' Calculating <nk|Sigma_c(omega)|nk> at k = ',kgw(:),ch10,&
 306 &  ' bands n = from ',ib1,' to ',ib2,ch10
 307  call wrtout(std_out,msg,'COLL')
 308 
 309  ABI_MALLOC(w_maxval,(minbnd:maxbnd))
 310  w_maxval = zero
 311 
 312  if ( ANY(gwc_ngfft(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,gwc_ngfft)
 313  gwc_mgfft   = MAXVAL(gwc_ngfft(1:3))
 314  gwc_fftalga = gwc_ngfft(7)/100 !; gwc_fftalgc=MOD(gwc_ngfft(7),10)
 315 
 316  if (Dtset%pawcross==1) mgfftf = MAXVAL(rho_ngfft(1:3))
 317 
 318  can_symmetrize = .FALSE.
 319  if (Sigp%symsigma>0) then
 320    can_symmetrize = .TRUE.
 321    if (Sigp%gwcalctyp >= 20) then
 322     do spin=1,Wfd%nsppol
 323       can_symmetrize(spin) = .not. esymm_failed(QP_sym(spin))
 324       if (.not.can_symmetrize(spin)) then
 325         write(msg,'(a,i0,4a)')&
 326 &         " Symmetrization cannot be performed for spin: ",spin,ch10,&
 327 &         " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg)
 328         MSG_WARNING(msg)
 329       end if
 330     end do
 331    end if
 332    if (wfd%nspinor == 2) MSG_WARNING("Symmetrization with nspinor=2 not implemented")
 333  end if
 334 
 335  ABI_UNUSED(Pawang%l_max)
 336 
 337  ! Print type of calculation.
 338  mod10=MOD(Sigp%gwcalctyp,10); sigma_type = sigma_type_from_key(mod10)
 339  call wrtout(std_out,sigma_type,'COLL')
 340 
 341  ! Set up logical flags for Sigma calculation
 342  if (mod10==SIG_GW_AC.and.Sigp%gwcalctyp/=1) then
 343    MSG_ERROR("not implemented")
 344  end if
 345 
 346  ! Initialize some values
 347  nspinor = Wfd%nspinor
 348  npwc = sigp%npwc
 349  spinor_padc(:,:)=RESHAPE([0, 0, npwc, npwc, 0, npwc, npwc, 0], [2, 4])
 350  ABI_MALLOC(npoles_missing,(minbnd:maxbnd))
 351  npoles_missing=0
 352 
 353  ! Normalization of theta_mu_minus_e0i
 354  ! If nsppol==2, qp_occ $\in [0,1]$
 355  SELECT CASE (Wfd%nsppol)
 356  CASE (1)
 357    fact_sp=half; tol_empty=0.01   ! below this value the state is assumed empty
 358    if (Wfd%nspinor==2) then
 359     fact_sp=one; tol_empty=0.005  ! below this value the state is assumed empty
 360    end if
 361  CASE (2)
 362    fact_sp=one; tol_empty=0.005   ! to be consistent and obtain similar results if a metallic
 363  CASE DEFAULT                     ! spin unpolarized system is treated using nsppol==2
 364    MSG_BUG('Wrong nsppol')
 365  END SELECT
 366 
 367  ! Allocate arrays used to accumulate the matrix elements of \Sigma_c over
 368  ! k-points and bands. Note that for AC requires only the imaginary frequencies
 369  !nomega_sigc=Sr%nomega_r+Sr%nomega4sd
 370  !if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i
 371  !
 372  ! === Define the G-G0 shifts for the FFT of the oscillators ===
 373  ! * Sigp%mG0 gives the MAX G0 component to account for umklapp.
 374  ! * Note the size MAX(Sigp%npwx,npwc).
 375  !
 376  ! === Precalculate the FFT index of $(R^{-1}(r-\tau))$ ===
 377  ! * S=\transpose R^{-1} and k_BZ = S k_IBZ
 378  ! * irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk.
 379  gwc_nfftot = PRODUCT(gwc_ngfft(1:3))
 380  ABI_MALLOC(irottb,(gwc_nfftot,Cryst%nsym))
 381  call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwc_ngfft,irottb,iscompatibleFFT)
 382  if (.not.iscompatibleFFT) then
 383    msg = "FFT mesh is not compatible with symmetries. Results might be affected by large errors!"
 384    MSG_WARNING(msg)
 385  end if
 386 
 387  ABI_MALLOC(ktabr,(gwc_nfftot,Kmesh%nbz))
 388  do ik_bz=1,Kmesh%nbz
 389    isym=Kmesh%tabo(ik_bz)
 390    do ifft=1,gwc_nfftot
 391      ktabr(ifft,ik_bz)=irottb(ifft,isym)
 392    end do
 393  end do
 394  ABI_FREE(irottb)
 395 
 396  if (Psps%usepaw==1 .and. Dtset%pawcross==1) then
 397    nfftf = PRODUCT(rho_ngfft(1:3))
 398    ABI_MALLOC(irottb,(nfftf,Cryst%nsym))
 399    call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,rho_ngfft,irottb,iscompatibleFFT)
 400 
 401    ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz))
 402    do ik_bz=1,Kmesh%nbz
 403      isym=Kmesh%tabo(ik_bz)
 404      do ifft=1,nfftf
 405        ktabrf(ifft,ik_bz)=irottb(ifft,isym)
 406      end do
 407    end do
 408    ABI_FREE(irottb)
 409  end if
 410 
 411  Sigcij_tab => Sigp%Sigcij_tab(ikcalc,1:Wfd%nsppol)
 412 
 413  got=0
 414  ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,Wfd%nsppol))
 415  call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,Wfd%nsppol,can_symmetrize,kgw,Sigp%mg0,&
 416 &  my_nbks,proc_distrb,got,global=.TRUE.)
 417 
 418  write(msg,'(a,i0,a)')" Will sum ",my_nbks," (b,k,s) states in Sigma_c."
 419  call wrtout(std_out,msg,'PERS')
 420 
 421  if (Sigp%gwcomp==1) then
 422    en_high=MAXVAL(qp_ene(Sigp%nbnds,:,:)) + Sigp%gwencomp
 423    write(msg,'(6a,e11.4,a)')ch10,&
 424 &    ' Using the extrapolar approximation to accelerate convergence',ch10,&
 425 &    ' with respect to the number of bands included',ch10,&
 426 &    ' with extrapolar energy: ',en_high*Ha_eV,' [eV]'
 427    call wrtout(std_out,msg,'COLL')
 428    ABI_MALLOC(wf1swf2_g,(gwc_nfftot*nspinor))
 429  endif
 430 
 431  if (Sigp%gwcomp == 1) then
 432    ! Setup of MPI table for extrapolar contributions.
 433    ABI_MALLOC(extrapolar_distrb,(ib1:ib2,ib1:ib2,Kmesh%nbz,Wfd%nsppol))
 434    extrapolar_distrb = xmpi_undefined_rank
 435 
 436    do spin=1,Wfd%nsppol
 437      do ik_bz=1,Kmesh%nbz
 438         if (ANY(proc_distrb(:,ik_bz,spin) /= xmpi_undefined_rank) ) then ! This BZ point will be calculated.
 439            rank_mask = .FALSE. ! The set of node that will treat (k,s).
 440            do band=1,Wfd%mband
 441              rank = proc_distrb(band,ik_bz,spin)
 442              if (rank /= xmpi_undefined_rank) rank_mask(rank+1)=.TRUE.
 443            end do
 444            do band2=ib1,ib2
 445              do irow=1,Sigcij_tab(spin)%col(band2)%size1   ! Looping over the non-zero elements of sigma_ij.
 446                band1 = Sigcij_tab(spin)%col(band2)%bidx(irow)
 447                idle = imin_loc(got,mask=rank_mask)
 448                got(idle) = got(idle)+1
 449                extrapolar_distrb(band1,band2,ik_bz,spin) = idle-1
 450              end do
 451            end do
 452         end if
 453      end do
 454    end do
 455 
 456    write(msg,'(a,i0,a)')" Will treat ",COUNT(extrapolar_distrb==Wfd%my_rank)," extrapolar terms."
 457    call wrtout(std_out,msg,'PERS')
 458  end if
 459 
 460  ABI_MALLOC(rhotwg_ki, (npwc*nspinor, minbnd:maxbnd))
 461  rhotwg_ki=czero_gw
 462  ABI_MALLOC(rhotwg, (npwc*nspinor))
 463  ABI_MALLOC(rhotwgp, (npwc*nspinor))
 464  ABI_MALLOC(vc_sqrt_qbz, (npwc))
 465 
 466  if (Er%mqmem==0) then ! Use out-of-core solution for epsilon.
 467    MSG_COMMENT('Reading q-slices from file. Slower but less memory.')
 468  end if                                                                                !
 469 
 470  ! Additional allocations for PAW.
 471  if (Psps%usepaw==1) then
 472    ABI_DT_MALLOC(Cprj_ksum,(Cryst%natom,nspinor))
 473    call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm)
 474    !
 475    ! For the extrapolar method we need the onsite terms of the PW in the FT mesh.
 476    ! * gw_gfft is the set of plane waves in the FFT Box for the oscillators.
 477    if (Sigp%gwcomp==1) then
 478      ABI_MALLOC(gw_gfft,(3,gwc_nfftot))
 479      q0=zero
 480      call get_gftt(gwc_ngfft,q0,Cryst%gmet,gw_gsq,gw_gfft)
 481      ABI_DT_MALLOC(Pwij_fft,(Psps%ntypat))
 482      call pawpwij_init(Pwij_fft,gwc_nfftot,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
 483    end if
 484  end if ! usepaw==1
 485 
 486  if (mod10==SIG_GW_AC) then ! Calculate Gauss-Legendre quadrature knots and weights for analytic continuation
 487    call coeffs_gausslegint(zero,one,gl_knots,gl_wts,Er%nomega_i)
 488 
 489    do io=1,Er%nomega_i ! First frequencies are always real
 490      if (ABS(AIMAG(one*Er%omega(Er%nomega_r+io))-(one/gl_knots(io)-one)) > 0.0001) then
 491       write(msg,'(3a)')&
 492 &       ' Frequencies in the SCR file are not compatible with the analytic continuation. ',ch10,&
 493 &       ' Verify the frequencies in the SCR file. '
 494       MSG_WARNING(msg)
 495       if (Wfd%my_rank==Wfd%master) write(std_out,*)AIMAG(Er%omega(Er%nomega_r+io)),(one/gl_knots(io)-one)
 496       MSG_ERROR("Cannot continue!")
 497      end if
 498    end do
 499 
 500    ! To calculate \int_0^\infty domegap f(omegap), we calculate \int_0^1 dz f(1/z-1)/z^2.
 501    omegap(:)=one/gl_knots(:)-one
 502    omegap2(:)=omegap(:)*omegap(:)
 503    ABI_MALLOC(ac_epsm1cqwz2, (npwc, npwc, Er%nomega_i))
 504    ABI_MALLOC(ac_integr, (npwc, npwc, Sr%nomega_i))
 505  end if
 506 
 507  ! Calculate total number of frequencies and allocate related arrays.
 508  ! sigcme2 is used to accumulate the diagonal matrix elements over k-points and
 509  ! GW bands, used only in case of ppmodel 3 and 4 (TODO save memory)
 510  nomega_tot=Sr%nomega_r+Sr%nomega4sd
 511  ABI_MALLOC(sigcme2,(nomega_tot,ib1:ib2))
 512  ABI_MALLOC(sigcme_3,(nomega_tot))
 513  sigcme2=czero_gw; sigcme_3=czero_gw
 514 
 515  ABI_MALLOC(sigctmp,(nomega_sigc,Sigp%nsig_ab))
 516  sigctmp=czero_gw
 517  ABI_MALLOC(sigc_ket, (npwc*nspinor, nomega_sigc))
 518 
 519  ! Arrays storing the contribution given by the Hermitian/anti-Hermitian part of \Sigma_c
 520  ABI_MALLOC(aherm_sigc_ket, (npwc*nspinor, nomega_sigc))
 521  ABI_MALLOC( herm_sigc_ket, (npwc*nspinor, nomega_sigc))
 522 
 523  sigcme_tmp=czero
 524 
 525  ABI_MALLOC(sigc,(2,nomega_sigc,ib1:ib2,ib1:ib2,Wfd%nsppol*Sigp%nsig_ab))
 526  sigc=czero
 527 
 528  !FIXME This quantities are only used for model GW if I am not wrong
 529  ABI_MALLOC(ket1, (npwc*nspinor, nomega_tot))
 530  ABI_MALLOC(ket2, (npwc*nspinor, nomega_tot))
 531  ABI_MALLOC(omegame0i,(nomega_tot))
 532 
 533  ! Here we divide the states where the QP energies are required into complexes. Note however that this approach is not
 534  ! based on group theory, and it might lead to spurious results in case of accidental degeneracies.
 535  nq_summed=Kmesh%nbz
 536  if (Sigp%symsigma>0) then
 537    call littlegroup_print(Ltg_k,std_out,Dtset%prtvol,'COLL')
 538    nq_summed=SUM(Ltg_k%ibzq(:))
 539    !
 540    ! Find number of degenerate subspaces and number of bands in each subspace
 541    ! The tolerance is a little bit arbitrary (0.001 eV)
 542    ! It could be reduced, in particular in case of nearly accidental degeneracies
 543    ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,Wfd%nsppol))
 544    degtab=0
 545    do spin=1,Wfd%nsppol
 546      do ib=ib1,ib2
 547        do jb=ib1,ib2
 548         if (ABS(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin))<0.001/Ha_ev) then
 549           degtab(ib,jb,spin)=1
 550         end if
 551        end do
 552      end do
 553    end do
 554  end if !symsigma
 555 
 556  write(msg,'(2a,i6,a)')ch10,' calculation status ( ',nq_summed,' to be completed):'
 557  call wrtout(std_out,msg,'COLL')
 558 
 559  ! Here we have a problem in case of CD since epsm1q might be huge
 560  ! TODO if single q (ex molecule) dont allocate epsm1q, avoid waste of memory
 561  if ( ANY(mod10== [SIG_GW_AC, SIG_GW_CD, SIG_QPGW_CD])) then
 562    if (.not.(mod10==SIG_GW_CD.and.Er%mqmem==0)) then
 563      ABI_STAT_MALLOC(epsm1_qbz, (npwc, npwc, Er%nomega), ierr)
 564      ABI_CHECK(ierr==0, "out-of-memory in epsm1_qbz")
 565    end if
 566  end if
 567 
 568  ! TODO epsm1_trcc_qbz is needed for SIG_GW_CD with symmetries since
 569  ! the Hermitian and the anti-Hermitian part have to be symmetrized in a different way.
 570  !if (mod10==SIG_QPGW_CD) then
 571  if (mod10==SIG_QPGW_CD) then
 572    ABI_STAT_MALLOC(epsm1_trcc_qbz, (npwc, npwc, Er%nomega), ierr)
 573    ABI_CHECK(ierr==0, "out-of-memory in epsm1_trcc_qbz")
 574    ABI_MALLOC(epsm1_tmp, (npwc, npwc))
 575  end if
 576 
 577  ABI_MALLOC(igfftcg0,(Gsph_Max%ng))
 578  ABI_MALLOC(ur_ibz,(gwc_nfftot*nspinor))
 579  ABI_MALLOC(usr_bz,(gwc_nfftot*nspinor))
 580 
 581  if (Dtset%pawcross==1) then
 582    ABI_MALLOC(igfftfcg0,(Gsph_c%ng))
 583    ABI_MALLOC(ur_ae_sum,(nfftf*nspinor))
 584    ABI_MALLOC(ur_ae_onsite_sum,(nfftf*nspinor))
 585    ABI_MALLOC(ur_ps_onsite_sum,(nfftf*nspinor))
 586  end if
 587  call timab(432,2,tsec) ! Init
 588 
 589  ! ==========================================
 590  ! ==== Fat loop over k_i in the full BZ ====
 591  ! ==========================================
 592  do spin=1,Wfd%nsppol
 593    if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE
 594    call timab(433,1,tsec) ! Init spin
 595 
 596    ! Load wavefunctions for GW corrections
 597    ! TODO: Rotate the functions here instead of calling rho_tw_g
 598    ABI_MALLOC(wfr_bdgw,(gwc_nfftot*nspinor,ib1:ib2))
 599    call wfd_get_many_ur(Wfd, [(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw)
 600 
 601    if (Wfd%usepaw==1) then
 602      ! Load cprj for GW states, note the indexing.
 603      dimcprj_gw=nspinor*(ib2-ib1+1)
 604      ABI_DT_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1))
 605      call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm)
 606      ibsp=ib1
 607      do jb=ib1,ib2
 608        call wfd_get_cprj(Wfd,jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
 609        call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
 610        call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1)))
 611        ibsp=ibsp+nspinor
 612      end do
 613      if (Dtset%pawcross==1) then
 614        ABI_MALLOC(ur_ae_bdgw,(nfftf*nspinor,ib1:ib2))
 615        ABI_MALLOC(ur_ae_onsite_bdgw,(nfftf*nspinor,ib1:ib2))
 616        ABI_MALLOC(ur_ps_onsite_bdgw,(nfftf*nspinor,ib1:ib2))
 617        do jb=ib1,ib2
 618          call wfd_paw_get_aeur(Wfdf,jb,jk_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,&
 619 &          ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
 620          ur_ae_bdgw(:,jb)=ur_ae_sum
 621          ur_ae_onsite_bdgw(:,jb)=ur_ae_onsite_sum
 622          ur_ps_onsite_bdgw(:,jb)=ur_ps_onsite_sum
 623        end do
 624      end if
 625    end if ! usepaw
 626 
 627    call timab(433,2,tsec) ! Init spin
 628 
 629    do ik_bz=1,Kmesh%nbz
 630      ! Parallelization over k-points and spin
 631      ! For the spin there is another check in the inner loop
 632      if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE
 633 
 634      call timab(434,1,tsec) ! initq
 635 
 636      ! Find the corresponding irreducible k-point
 637      call get_BZ_item(Kmesh,ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt)
 638      spinrot_kbz(:)=Cryst%spinrot(:,isym_ki)
 639 
 640      ! Identify q and G0 where q + G0 = k_GW - k_i
 641      kgw_m_ksum=kgw-ksum
 642      call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0)
 643 
 644      ! If symsigma, symmetrize the matrix elements.
 645      ! Sum only q"s in IBZ_k. In this case elements are weighted
 646      ! according to wtqp and wtqm. wtqm is for time-reversal.
 647      wtqp=1; wtqm=0
 648      if (can_symmetrize(spin)) then
 649        if (Ltg_k%ibzq(iq_bz)/=1) CYCLE
 650        wtqp=0; wtqm=0
 651        do isym=1,Ltg_k%nsym_sg
 652          wtqp=wtqp+Ltg_k%wtksym(1,isym,iq_bz)
 653          wtqm=wtqm+Ltg_k%wtksym(2,isym,iq_bz)
 654        end do
 655      end if
 656 
 657      write(msg,'(3(a,i0),a,i0)')' Sigma_c: ik_bz ',ik_bz,'/',Kmesh%nbz,", spin: ",spin,' done by mpi-rank: ',Wfd%my_rank
 658      call wrtout(std_out,msg,'PERS')
 659 
 660      ! Find the corresponding irred q-point.
 661      call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q)
 662      q_is_gamma = (normv(qbz,Cryst%gmet,"G") < GW_TOL_W0)
 663 
 664      !q_is_gamma = (normv(qbz,Cryst%gmet,"G") < 0.7)
 665      !if (iq_ibz/=2.and.iq_ibz/=1) CYCLE
 666      !if (ANY(qbz<=-(half-tol16)) .or. ANY(qbz>(half+tol16))) CYCLE
 667      !if (q_is_gamma) then; write(std_out,*)"skipping q=Gamma"; CYCLE; end if
 668      !
 669      ! Tables for the FFT of the oscillators.
 670      !  a) FFT index of the G-G0.
 671      !  b) gw_gbound table for the zero-padded FFT performed in rhotwg.
 672      ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2))
 673      call gsph_fft_tabs(Gsph_c,g0,gwc_mgfft,gwc_ngfft,use_padfft,gw_gbound,igfftcg0)
 674 
 675      if (ANY(gwc_fftalga == [2, 4])) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
 676 #ifdef FC_IBM
 677      ! XLF does not deserve this optimization (problem with [v67mbpt][t03])
 678      use_padfft = 0
 679 #endif
 680      if (use_padfft==0) then
 681        ABI_FREE(gw_gbound)
 682        ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2*use_padfft))
 683      end if
 684 
 685      if (Dtset%pawcross==1) then
 686        ABI_MALLOC(gboundf,(2*mgfftf+8,2))
 687        call gsph_fft_tabs(Gsph_c,g0,mgfftf,rho_ngfft,use_padfftf,gboundf,igfftfcg0)
 688        if ( ANY(gwc_fftalga == (/2,4/)) ) use_padfftf=0
 689        if (use_padfftf==0) then
 690          ABI_FREE(gboundf)
 691          ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf))
 692        end if
 693      end if
 694 
 695      ! Evaluate oscillator matrix elements
 696      ! $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form.
 697      if (Psps%usepaw==1) then
 698        ABI_DT_MALLOC(Pwij_qg,(Psps%ntypat))
 699        q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT
 700        call pawpwij_init(Pwij_qg,npwc,q0,Gsph_c%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
 701      end if
 702 
 703      if (Er%mqmem==0) then
 704        ! Read q-slice of epsilon^{-1}|chi0 in Er%epsm1(:,:,:,1) (much slower but less memory).
 705        call get_epsm1(Er,Vcp,0,0,Dtset%iomode,xmpi_comm_self,iqibzA=iq_ibz)
 706        if (sigma_needs_ppm(Sigp)) then
 707          if (Wfd%usepaw==1.and.PPm%userho==1) then
 708            ! Use PAW AE rhor.
 709            call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,rho_nfftot,Gsph_c%gvec,&
 710 &            rho_ngfft,aepaw_rhor(:,1),iqiA=iq_ibz)
 711          else
 712            call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,rho_nfftot,Gsph_c%gvec,&
 713 &            rho_ngfft,rhor(:,1),iqiA=iq_ibz)
 714          end if
 715        end if
 716      end if
 717 
 718      ! Symmetrize PPM parameters and epsm1 (q_IBZ --> q_BZ):
 719      ! NOTE:
 720      !    - We are not considering umklapp with G0/=0. In this case,
 721      !      indeed the equation is different since we have to use G-G0.
 722      !      A check, however, is performed in sigma.
 723      !    - If gwcomp==1 and mod10 in [1,2,9], one needs both to set up botsq and epsm1_q
 724      if (sigma_needs_ppm(Sigp)) then
 725        ABI_MALLOC(botsq,(PPm%npwc,PPm%dm2_botsq))
 726        ABI_MALLOC(otq,(PPm%npwc,PPm%dm2_otq))
 727        ABI_MALLOC(eig,(PPm%dm_eig,PPm%dm_eig))
 728        call ppm_get_qbz(PPm,Gsph_c,Qmesh,iq_bz,botsq,otq,eig)
 729      end if
 730 
 731      if ( ANY(mod10 == [SIG_GW_AC,SIG_GW_CD,SIG_QPGW_CD])) then
 732        ! Numerical integration or model GW with contour deformation or Analytic Continuation
 733        ! TODO In case of AC we should symmetrize only the imaginary frequencies
 734        if (mod10==SIG_GW_CD.and.Er%mqmem==0) then
 735          ! Do in-place symmetrisation
 736          call Epsm1_symmetrizer_inplace(iq_bz,Er%nomega, npwc,Er,Gsph_c,Qmesh,.TRUE.)
 737        else
 738          call Epsm1_symmetrizer(iq_bz,Er%nomega, npwc,Er,Gsph_c,Qmesh,.TRUE.,epsm1_qbz)
 739        end if
 740 
 741        if (mod10==SIG_GW_AC) then
 742          ! Prepare the first term: \Sum w_i 1/z_i^2 f(1/z_i-1)..
 743          ! The first frequencies are always real, skip them.
 744          ! Memory is not optimized.
 745          do iiw=1,Er%nomega_i
 746            z2=gl_knots(iiw)*gl_knots(iiw)
 747            ac_epsm1cqwz2(:,:,iiw)= gl_wts(iiw)*epsm1_qbz(:,:,Er%nomega_r+iiw)/z2
 748          end do
 749        end if
 750 
 751        if (mod10==SIG_QPGW_CD) then
 752          ! For model GW we need transpose(conjg(epsm1_qbz)) ===
 753          do io=1,Er%nomega
 754           epsm1_tmp(:,:) = GWPC_CONJG(epsm1_qbz(:,:,io))
 755           epsm1_trcc_qbz(:,:,io) = TRANSPOSE(epsm1_tmp)
 756          end do
 757        end if
 758      end if !gwcalctyp
 759 
 760      ! Get Fourier components of the Coulombian interaction in the BZ ===
 761      ! In 3D systems, neglecting umklapp: vc(Sq,sG) = vc(q,G) = 4pi/|q+G|**2
 762      ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
 763      do ig=1,npwc
 764        vc_sqrt_qbz(Gsph_c%rottb(ig,itim_q,isym_q))=Vcp%vc_sqrt(ig,iq_ibz)
 765      end do
 766 
 767      call timab(434,2,tsec) ! initq
 768 
 769      ! Sum over band
 770      call timab(445,1,tsec) ! loop
 771      do ib=1,Sigp%nbnds
 772        call timab(436,1,tsec) ! (1)
 773 
 774        ! Parallelism over spin
 775        ! This processor has this k-point but what about spin?
 776        if (proc_distrb(ib,ik_bz,spin)/=Wfd%my_rank) CYCLE
 777 
 778        call wfd_get_ur(Wfd,ib,ik_ibz,spin,ur_ibz)
 779 
 780        if (Psps%usepaw==1) then
 781          ! Load cprj for point ksum, this spin or spinor and *THIS* band.
 782          ! TODO MG I could avoid doing this but I have to exchange spin and bands ???
 783          ! For sure there is a better way to do this!
 784          call wfd_get_cprj(Wfd,ib,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
 785          call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
 786          if (Dtset%pawcross==1) then
 787            call wfd_paw_get_aeur(Wfdf,ib,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,&
 788 &              ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
 789          end if
 790        end if
 791 
 792        if (mod10==SIG_GW_AC) then
 793          ! Calculate integral over omegap with Gauss-Legendre quadrature.
 794          ! * -1/pi \int domegap epsm1c*(omega-e0i) / ( (omega-e0i)^2 + omegap^2)
 795          ! * Note that energies are calculated wrt the Fermi level.
 796          ac_integr(:,:,:)=czero_gw
 797          do io=1,Sr%nomega_i
 798            omegame0i_ac  = Sr%omega_i(io)-qp_ene(ib,ik_ibz,spin)
 799            omegame0i2_ac = omegame0i_ac*omegame0i_ac
 800            do iiw=1,Er%nomega_i
 801              do iggp=0,npwc*npwc-1
 802                ig=iggp/npwc+1
 803                igp= iggp-(ig-1)*npwc+1 ! \int domegap epsm1c/((omega-e0i)^2 + omegap^2)
 804                ac_integr(ig,igp,io)= ac_integr(ig,igp,io) + ac_epsm1cqwz2(ig,igp,iiw)/(omegame0i2_ac + omegap2(iiw))
 805              end do
 806            end do
 807            ac_integr(:,:,io)=ac_integr(:,:,io)*omegame0i_ac
 808          end do
 809          ac_integr(:,:,:)=-ac_integr(:,:,:)*piinv
 810        end if
 811 
 812        call timab(436,2,tsec) ! (1)
 813        call timab(437,1,tsec) ! rho_tw_g
 814 
 815        ! Get all <k-q,ib,s|e^{-i(q+G).r}|s,jb,k>, at once
 816        do jb=ib1,ib2
 817 
 818          call rho_tw_g(nspinor,npwc,gwc_nfftot,ndat1,gwc_ngfft,1,use_padfft,igfftcg0,gw_gbound,&
 819 &          ur_ibz        ,iik,ktabr(:,ik_bz),ph_mkt  ,spinrot_kbz,  &
 820 &          wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,&
 821 &          nspinor,rhotwg_ki(:,jb))
 822 
 823          if (Psps%usepaw==1) then
 824            ! Add on-site contribution, projectors are already in BZ !TODO Recheck this!
 825            i2=jb; if (nspinor==2) i2=(2*jb-1)
 826            spad=(nspinor-1)
 827            call paw_rho_tw_g(npwc,nspinor,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_c%gvec,&
 828 &            Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),Pwij_qg,rhotwg_ki(:,jb))
 829            if (Dtset%pawcross==1) then ! Add paw cross term
 830              call paw_cross_rho_tw_g(nspinor,npwc,nfftf,rho_ngfft,1,use_padfftf,igfftfcg0,gboundf,&
 831 &             ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum,iik,ktabrf(:,ik_bz),ph_mkt,spinrot_kbz,&
 832 &             ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,&
 833 &             nspinor,rhotwg_ki(:,jb))
 834            end if
 835          end if
 836 
 837          ! Multiply by the square root of the Coulomb term
 838          ! In 3-D systems, the factor sqrt(4pi) is included)
 839          do ii=1,nspinor
 840            spad=(ii-1)*npwc
 841            rhotwg_ki(spad+1:spad+npwc,jb) = rhotwg_ki(spad+1:spad+npwc,jb)*vc_sqrt_qbz(1:npwc)
 842          end do
 843 
 844          ! === Treat analytically the case q --> 0 ===
 845          ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma
 846          !   while the Colulomb term is integrated out.
 847          ! * In the scalar case we have nonzero contribution only if ib==jb
 848          ! * For nspinor==2 evalute <ib,up|jb,up> and <ib,dwn|jb,dwn>,
 849          !   impose orthonormalization since npwwfn might be < npwvec.
 850          if (ik_bz==jk_bz) then
 851            if (nspinor==1) then
 852              rhotwg_ki(1,jb)=czero_gw
 853              if (ib==jb) rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp)
 854 
 855            else
 856              npw_k = Wfd%npwarr(ik_ibz)
 857              rhotwg_ki(1, jb) = zero; rhotwg_ki(npwc+1, jb) = zero
 858              if (ib==jb) then
 859                cg_sum => Wfd%Wave(ib,ik_ibz,spin)%ug
 860                cg_jb  => Wfd%Wave(jb,jk_ibz,spin)%ug
 861                ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1)
 862                rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
 863                ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1)
 864                rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
 865                ! PAW is missing
 866              end if
 867            end if
 868          end if
 869        end do !jb  Got all matrix elements from minbnd up to maxbnd.
 870 
 871        theta_mu_minus_e0i=fact_sp*qp_occ(ib,ik_ibz,spin)
 872 
 873        ! Starting point to evaluate the derivative of Sigma and the Spectral function
 874        e0i=qp_ene(ib,ik_ibz,spin)
 875 
 876        ! Frequencies for the spectral function, e0i=qp_ene(ib,ik_ibz,spin)
 877        ! FIXME the interval is not centered on eoi ! WHY?
 878        if (Sr%nomega_r>0) omegame0i(1:Sr%nomega_r)=DBLE(Sr%omega_r(1:Sr%nomega_r))-e0i
 879 
 880        call timab(437,2,tsec) ! rho_tw_g
 881 
 882        do kb=ib1,ib2
 883          call timab(438,1,tsec) ! (2)
 884 
 885          ! Get frequencies $\omega$-\epsilon_in$ to evaluate  $d\Sigma/dE$, note the spin
 886          ! subtract e_KS since we have stored e_KS+ Delta \omega in Sr%omega4sd, not required for AC
 887          do io=Sr%nomega_r+1,nomega_tot
 888            omegame0i(io)=DBLE(Sr%omega4sd(kb,jk_ibz,io-Sr%nomega_r,spin))-e0i
 889          end do
 890 
 891          ! Get the ket \Sigma|\phi_{k,kb}> according to the method.
 892          rhotwgp(:)=rhotwg_ki(:,kb)
 893          sigc_ket  = czero_gw
 894          ket1      = czero_gw
 895          ket2      = czero_gw
 896          !aherm_sigc_ket = czero_gw
 897          ! herm_sigc_ket = czero_gw
 898 
 899          SELECT CASE (mod10)
 900          CASE (SIG_GW_PPM)
 901            ! GW WITH Plasmon-Pole Model.
 902            ! Note that ppmodel 3 or 4 work only in case of standard perturbative approach!
 903            ! Moreover, for ppmodel 3 and 4, spinorial case is not allowed
 904            call calc_sig_ppm(PPm,nspinor,npwc,nomega_tot,rhotwgp,botsq,otq,&
 905 &           omegame0i,Sigp%zcut,theta_mu_minus_e0i,eig,npwc,sigc_ket,sigcme_3)
 906 
 907            if (PPm%model==3.or.PPm%model==4) then
 908              sigcme2(:,kb)=sigcme2(:,kb) + (wtqp+wtqm)*DBLE(sigcme_3(:)) + (wtqp-wtqm)*j_gw*AIMAG(sigcme_3(:))
 909            end if
 910 
 911          CASE (SIG_GW_AC)
 912            ! GW with Analytic continuation.
 913            ! Evaluate \sum_Gp integr_GGp(omegasi) rhotw_Gp TODO this part can be optimized
 914            do io=1,Sr%nomega_i
 915              do ispinor=1,nspinor
 916                spadc=(ispinor-1)*npwc
 917                do ig=1,npwc
 918                  ctmp=czero
 919                  do igp=1,npwc
 920                    ctmp=ctmp+ac_integr(ig,igp,io)*rhotwgp(igp+spadc)
 921                  end do
 922                  sigc_ket(ig+spadc,io)=ctmp
 923                end do
 924              end do !ispinor
 925            end do !io
 926 
 927          CASE (SIG_GW_CD)
 928            ! GW with contour deformation.
 929            ! Check if pole contributions need to be summed
 930            ! (this avoids unnecessary splint calls and saves time)
 931            !me_calc_poles = .TRUE.
 932            do io=1,nomega_tot
 933              if (omegame0i(io)>=zero.AND.(ABS(one-theta_mu_minus_e0i)>zero)) then
 934                !me_calc_poles(io) = .TRUE.
 935                if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io))
 936              else if (omegame0i(io)<zero.AND.(ABS(theta_mu_minus_e0i)>zero)) then
 937                !me_calc_poles(io) = .TRUE.
 938                if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io))
 939              end if
 940            end do
 941 
 942            ! Check memory saving
 943            if (Er%mqmem == 0) then
 944              call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,&
 945 &             Er%omega,Er%epsm1(:,:,:,1),omegame0i,theta_mu_minus_e0i,sigc_ket,Dtset%ppmfrq,npoles_missing(kb),&
 946 &              method=Dtset%cd_frqim_method)
 947            else
 948              call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,&
 949 &             Er%omega,epsm1_qbz,omegame0i,theta_mu_minus_e0i,sigc_ket,Dtset%ppmfrq,npoles_missing(kb),&
 950 &              method=Dtset%cd_frqim_method)
 951            end if
 952 
 953 #if 0
 954            if (wtqm/=0) then
 955              call calc_sigc_cd(npwc,npwc,nspinor,,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,&
 956 &              Er%omega,epsm1_trcc_qbz,omegame0i,theta_mu_minus_e0i,aherm_sigc_ket,Dtset%ppmfrq,npoles_missing(kb),&
 957 &               method=Dtset%cd_frqim_method)
 958 
 959               herm_sigc_ket = half*(sigc_ket + aherm_sigc_ket)
 960              aherm_sigc_ket = half*(sigc_ket - aherm_sigc_ket)
 961            else
 962              herm_sigc_ket  = sigc_ket
 963              aherm_sigc_ket = czero_gw
 964            end if
 965 #endif
 966 
 967          CASE (SIG_QPGW_PPM)
 968            ! MODEL GW calculation WITH PPm  TODO Spinor not tested.
 969            ! Calculate \Sigma(E_k) |k> to obtain <j|\Sigma(E_k)|k>
 970            ABI_MALLOC(sigcme_new,(nomega_tot))
 971 
 972            call calc_sig_ppm(PPm,nspinor,npwc,nomega_tot,rhotwgp,botsq,otq,&
 973 &            omegame0i,Sigp%zcut,theta_mu_minus_e0i,eig,npwc,ket1,sigcme_new)
 974 
 975            if (Sigp%gwcalctyp==28) then
 976              if (PPm%model/=1.and.PPm%model/=2) then
 977                ! This is needed to have npwc=PPm%dm2_botsq=PPm%dm2_otq
 978                write(msg,'(3a)')&
 979 &                'For the time being, gwcalctyp=28 cannot be used with ppmodel=3,4.',ch10,&
 980 &                'Use another Plasmon Pole Model when gwcalctyp=28.'
 981                MSG_ERROR(msg)
 982              end if
 983              ABI_MALLOC(botsq_conjg_transp,(PPm%dm2_botsq,npwc))
 984              botsq_conjg_transp=TRANSPOSE(botsq) ! Keep these two lines separated, otherwise gfortran messes up
 985              botsq_conjg_transp=CONJG(botsq_conjg_transp)
 986              ABI_MALLOC(otq_transp,(PPm%dm2_otq,PPm%npwc))
 987              otq_transp=TRANSPOSE(otq)
 988 
 989              call calc_sig_ppm(PPm,nspinor,npwc,nomega_tot,rhotwgp,botsq_conjg_transp,otq_transp,&
 990 &              omegame0i,Sigp%zcut,theta_mu_minus_e0i,eig,npwc,ket2,sigcme_3)
 991 
 992              ABI_FREE(botsq_conjg_transp)
 993              ABI_FREE(otq_transp)
 994              sigc_ket= half*(ket1+ket2)
 995            else
 996              sigc_ket= ket1
 997            end if
 998 
 999            ABI_FREE(sigcme_new)
1000 
1001          CASE (SIG_QPGW_CD)
1002            ! MODEL GW with numerical integration.
1003            ! Check if pole contributions need to be summed
1004            ! (this avoids unnecessary splint calls and saves time)
1005            !me_calc_poles = .TRUE.
1006            do io=1,nomega_tot
1007              if (omegame0i(io)>=zero.AND.(ABS(one-theta_mu_minus_e0i)>zero)) then
1008                !me_calc_poles(io) = .TRUE.
1009                if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io))
1010              else if (omegame0i(io)<zero.AND.(ABS(theta_mu_minus_e0i)>zero)) then
1011                !me_calc_poles(io) = .TRUE.
1012                if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io))
1013              end if
1014            end do
1015 
1016            ! Calculate \Sigma(E_k)|k> to obtain <j|\Sigma(E_k)|k>
1017            call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,&
1018 &            Er%omega,epsm1_qbz,omegame0i,theta_mu_minus_e0i,ket1,Dtset%ppmfrq,npoles_missing(kb),&
1019 &             method=Dtset%cd_frqim_method)
1020 
1021            if (Sigp%gwcalctyp==29) then
1022              ! Calculate \Sigma^*(E_k)|k> to obtain <k|\Sigma(E_k)|j>^*
1023              call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,&
1024 &              Er%omega,epsm1_trcc_qbz,omegame0i,theta_mu_minus_e0i,ket2,Dtset%ppmfrq,npoles_missing(kb),&
1025 &               method=Dtset%cd_frqim_method)
1026              sigc_ket = half*(ket1+ket2)
1027            else
1028              sigc_ket = ket1
1029            end if
1030 
1031          CASE DEFAULT
1032            MSG_ERROR(sjoin("Unsupported value for mod10:", itoa(mod10)))
1033          END SELECT
1034 
1035          if (Sigp%gwcomp==1) then
1036            ! TODO spinor not implemented
1037            call calc_sig_ppm_comp(npwc,nomega_tot,rhotwgp,botsq,otq,DBLE(Sr%egw(kb,jk_ibz,spin)-en_high),&
1038 &            Sigp%zcut,theta_mu_minus_e0i,sigc_ket,PPm%model,npwc,PPm%dm2_botsq,PPm%dm2_otq)
1039          end if
1040 
1041          call timab(438,2,tsec) !
1042          call timab(439,1,tsec) ! sigma_me
1043 
1044          ! Loop over the non-zero row elements of this column.
1045          ! 1) If gwcalctyp<20 : only diagonal elements since QP==KS.
1046          ! 2) If gwcalctyp>=20: only off-diagonal elements connecting states with same character.
1047          do irow=1,Sigcij_tab(spin)%col(kb)%size1
1048            jb = Sigcij_tab(spin)%col(kb)%bidx(irow)
1049            rhotwg=rhotwg_ki(:,jb)
1050 
1051            ! Calculate <\phi_j|\Sigma_c|\phi_k>
1052            ! Different freqs according to method (AC or Perturbative), see nomega_sigc.
1053            do iab=1,Sigp%nsig_ab
1054              spadc1 = spinor_padc(1, iab); spadc2 = spinor_padc(2, iab)
1055              do io=1,nomega_sigc
1056                sigctmp(io,iab) = XDOTC(npwc,rhotwg(spadc1+1:),1,sigc_ket(spadc2+1:,io),1)
1057              end do
1058            end do
1059 
1060            if (Sigp%gwcomp==1) then
1061              ! Evaluate Extrapolar term TODO this does not work with spinor
1062              if (extrapolar_distrb(jb,kb,ik_bz,spin) == Wfd%my_rank ) then
1063                ! Do it once as it does not depend on the ib index summed over.
1064                extrapolar_distrb(jb,kb,ik_bz,spin) = xmpi_undefined_rank
1065 #if 1
1066                call calc_wfwfg(ktabr(:,jk_ibz),jik, spinrot_kgw, & ! TODO: why jk_ibz?
1067                  gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g)
1068 #else
1069                call calc_wfwfg(ktabr(:,jk_bz),jik,spinrot_kgw,&
1070                  gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g)
1071 #endif
1072 
1073                if (Psps%usepaw==1) then
1074                  i1=jb; i2=kb
1075                  if (nspinor==2) then
1076                    i1=(2*jb-1); i2=(2*kb-1)
1077                  end if
1078                  spad=(nspinor-1)
1079                  call paw_rho_tw_g(gwc_nfftot,Sigp%nsig_ab,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,&
1080 &                  gw_gfft,Cprj_kgw(:,i1:i1+spad),Cprj_kgw(:,i2:i2+spad),Pwij_fft,wf1swf2_g)
1081 
1082                  if (Dtset%pawcross==1) then ! Add paw cross term
1083                    call paw_cross_rho_tw_g(nspinor,npwc,nfftf,rho_ngfft,1,use_padfftf,igfftfcg0,gboundf,&
1084 &                   ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,&
1085 &                   ur_ae_bdgw(:,kb),ur_ae_onsite_bdgw(:,kb),ur_ps_onsite_bdgw(:,kb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,&
1086 &                   nspinor,wf1swf2_g)
1087                  end if
1088                end if
1089 
1090                ! The static contribution from completeness relation is calculated once ===
1091                call calc_coh_comp(iq_ibz,Vcp%i_sz,(jb==kb),nspinor,Sigp%nsig_ab,DBLE(Sr%egw(kb,jk_ibz,spin)-en_high),&
1092 &                npwc,Gsph_c%gvec,gwc_ngfft,gwc_nfftot,wf1swf2_g,vc_sqrt_qbz,botsq,otq,sigcohme)
1093 
1094                do io=1,nomega_sigc
1095                  sigctmp(io,:) = sigctmp(io,:)+sigcohme(:)
1096                end do
1097              end if ! gwcomp==1
1098            end if ! gwcom==1
1099 
1100            ! Accumulate and, in case, symmetrize matrix elements of Sigma_c
1101            do iab=1,Sigp%nsig_ab
1102              is_idx=spin; if (nspinor==2) is_idx=iab
1103 
1104              sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + &
1105 &              (wtqp+wtqm)*DBLE(sigctmp(:,iab)) + (wtqp-wtqm)*j_gw*AIMAG(sigctmp(:,iab))
1106 
1107              sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + wtqp*      sigctmp(:,iab)
1108              sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + wtqm*CONJG(sigctmp(:,iab))
1109              ! TODO this should be the contribution coming from the anti-hermitian part.
1110            end do
1111          end do ! irow used to calculate matrix elements of $\Sigma$
1112 
1113          ! shaltaf (030406): this has to be done in a clean way later.
1114          ! TODO does not work with spinor.
1115          if (mod10==SIG_GW_PPM.and.(PPm%model==3.or.PPm%model==4)) then
1116            sigcme_tmp(:,kb,kb,spin)= sigcme2(:,kb)
1117            sigc(1,:,kb,kb,spin)= sigcme2(:,kb)
1118            sigc(2,:,kb,kb,spin)= czero
1119          end if
1120 
1121          call timab(439,2,tsec) ! csigme(SigC)
1122        end do !kb to calculate matrix elements of $\Sigma$
1123      end do !ib
1124 
1125      call timab(445,2,tsec) ! csigme(SigC)
1126 
1127      ! Deallocate k-dependent quantities.
1128      ABI_FREE(gw_gbound)
1129      if (Dtset%pawcross==1) then
1130        ABI_FREE(gboundf)
1131      end if
1132 
1133      if (sigma_needs_ppm(Sigp)) then
1134        ABI_FREE(botsq)
1135        ABI_FREE(otq)
1136        ABI_FREE(eig)
1137      end if
1138      if (Psps%usepaw==1) then
1139        call pawpwij_free(Pwij_qg)
1140        ABI_DT_FREE(Pwij_qg)
1141      end if
1142 
1143    end do ! ik_bz
1144 
1145    ABI_FREE(wfr_bdgw)
1146    if (Wfd%usepaw==1) then
1147      call pawcprj_free(Cprj_kgw)
1148      ABI_DT_FREE(Cprj_kgw)
1149      if (Dtset%pawcross==1) then
1150        ABI_FREE(ur_ae_bdgw)
1151        ABI_FREE(ur_ae_onsite_bdgw)
1152        ABI_FREE(ur_ps_onsite_bdgw)
1153      end if
1154    end if
1155  end do ! spin
1156 
1157  ABI_FREE(sigcme2)
1158  ABI_FREE(sigcme_3)
1159  ABI_FREE(igfftcg0)
1160  if (Dtset%pawcross==1) then
1161    ABI_FREE(igfftfcg0)
1162  end if
1163 
1164  ! Gather contributions from all the CPUs
1165  call timab(440,1,tsec) ! wfd_barrier
1166  call timab(440,2,tsec) ! wfd_barrier
1167  call timab(441,1,tsec) ! xmpi_sum
1168 
1169  call xmpi_sum(sigcme_tmp, wfd%comm, ierr)
1170  call xmpi_sum(sigc, wfd%comm, ierr)
1171  call timab(441,2,tsec) ! xmpi_sum
1172 
1173  ! Multiply by constants. In 3D systems sqrt(4pi) is included in vc_sqrt_qbz ===
1174  sigcme_tmp = sigcme_tmp /(Cryst%ucvol*Kmesh%nbz)
1175  sigc       = sigc       /(Cryst%ucvol*Kmesh%nbz)
1176 
1177  ! If we have summed over the IBZ_q now we have to average over complexes ===
1178  ! Presently only diagonal terms are considered
1179  ! TODO QP-SCGW required a more involved approach, there is a check in sigma
1180  ! TODO it does not work if nspinor==2.
1181  call timab(442,1,tsec) ! final ops
1182 
1183  do spin=1,Wfd%nsppol
1184    if (can_symmetrize(spin)) then
1185      if (mod10==SIG_GW_AC) then ! FIXME here there is a problem in case of AC with symmetries
1186        ABI_MALLOC(sym_cme, (Sr%nomega_i, ib1:ib2, ib1:ib2, sigp%nsig_ab))
1187      else
1188        ABI_MALLOC(sym_cme, (nomega_tot, ib1:ib2, ib1:ib2, sigp%nsig_ab))
1189      end if
1190      sym_cme=czero
1191 
1192      ! Average over degenerate diagonal elements
1193      ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results.
1194      ! another good reason to use a strict criterion for the tolerance on eigenvalues.
1195      do ib=ib1,ib2
1196        ndegs=0
1197        do jb=ib1,ib2
1198          if (degtab(ib,jb,spin)==1) then
1199            if (nspinor == 1) then
1200              sym_cme(:, ib, ib, 1) = sym_cme(:, ib, ib, 1) + SUM(sigc(:,:,jb,jb,spin), DIM=1)
1201            else
1202              do ii=1,sigp%nsig_ab
1203                sym_cme(:, ib, ib, ii) = sym_cme(:, ib, ib, ii) + SUM(sigc(:,:,jb,jb,ii), dim=1)
1204              end do
1205            end if
1206          end if
1207          ndegs = ndegs + degtab(ib,jb,spin)
1208        end do
1209        sym_cme(:,ib,ib,:) = sym_cme(:,ib,ib,:) / ndegs
1210      end do
1211 
1212      if (Sigp%gwcalctyp >= 20) then
1213        do iwc=1,nomega_sigc
1214          call esymm_symmetrize_mels(QP_sym(spin),ib1,ib2,sigc(:,iwc,:,:,spin),sym_cme(iwc,:,:,1))
1215        end do
1216      end if
1217 
1218      ! Copy symmetrized values
1219      do ib=ib1,ib2
1220        do jb=ib1,ib2
1221          !if (mod10==SIG_GW_AC.and.average_real) CYCLE ! this is to check another scheme in case of AC
1222          if (nspinor == 1) then
1223            sigcme_tmp(:,ib,jb,spin) = sym_cme(:,ib,jb,1)
1224          else
1225            do ii=1,sigp%nsig_ab
1226              sigcme_tmp(:,ib,jb,ii) = sym_cme(:,ib,jb,ii)
1227            end do
1228          end if
1229        end do
1230      end do
1231      ABI_FREE(sym_cme)
1232    end if
1233  end do
1234 
1235  ! Reconstruct the full sigma matrix from the upper triangle (only for HF, SEX and COHSEX)
1236  !if (Sigp%gwcalctyp>=20 .and. sigma_is_herm(Sigp) ) then
1237  !  ABI_CHECK(nspinor==1,"cannot hermitianize non-collinear sigma!")
1238  !  do spin=1,Wfd%nsppol
1239  !    do io=1,nomega_sigc
1240  !      call hermitianize(sigcme_tmp(io,:,:,spin),"Upper")
1241  !    end do
1242  !  end do
1243  !end if
1244 
1245  ! GW with contour deformation: check on the number of poles not included.
1246  if (ANY(mod10 == [SIG_GW_CD,SIG_QPGW_CD])) then
1247    call xmpi_sum(npoles_missing, wfd%comm, ierr)
1248    npls = SUM(npoles_missing)
1249    if (npls>0) then
1250      MSG_WARNING(sjoin("Total number of missing poles for contour deformation method:", itoa(npls)))
1251      do band=minbnd,maxbnd
1252        npls = npoles_missing(band)
1253        if (npls>0) then
1254          write(msg,'(a,2(i0,a))')" For band ",band," there are ",npls," missing poles"
1255          call wrtout(std_out,msg,"COLL")
1256        end if
1257      end do
1258    end if
1259    ! Print data on the maximum value needed for the screening along the real axis
1260    w_localmax = MAXVAL(w_maxval)
1261    call xmpi_max(w_localmax,w_max, wfd%comm, ierr)
1262    write(msg,'(a,f12.5,a)') ' Max omega value used in W(omega): ',w_max*Ha_eV,' [eV]'
1263    call wrtout(std_out,msg,"COLL")
1264  end if
1265  call timab(442,2,tsec) ! final ops
1266 
1267  ! ===========================
1268  ! ==== Deallocate memory ====
1269  ! ===========================
1270  if (Psps%usepaw==1) then
1271    if (allocated(gw_gfft)) then
1272      ABI_FREE(gw_gfft)
1273    end if
1274    call pawcprj_free(Cprj_ksum)
1275    ABI_DT_FREE(Cprj_ksum)
1276    if (allocated(Pwij_fft)) then
1277      call pawpwij_free(Pwij_fft)
1278      ABI_DT_FREE(Pwij_fft)
1279    end if
1280    if (Dtset%pawcross==1) then
1281      ABI_FREE(ur_ae_sum)
1282      ABI_FREE(ur_ae_onsite_sum)
1283      ABI_FREE(ur_ps_onsite_sum)
1284      ABI_FREE(ktabrf)
1285    end if
1286  end if
1287 
1288  ABI_FREE(npoles_missing)
1289  ABI_FREE(ur_ibz)
1290  ABI_FREE(usr_bz)
1291  ABI_FREE(ktabr)
1292  ABI_FREE(sigc_ket)
1293  ABI_FREE(ket1)
1294  ABI_FREE(ket2)
1295  ABI_FREE(rhotwg_ki)
1296  ABI_FREE(rhotwg)
1297  ABI_FREE(rhotwgp)
1298  ABI_FREE(vc_sqrt_qbz)
1299  ABI_FREE(omegame0i)
1300  ABI_FREE(sigctmp)
1301  ABI_FREE(sigc)
1302  ABI_FREE(w_maxval)
1303 
1304  if (allocated(epsm1_qbz)) then
1305    ABI_FREE(epsm1_qbz)
1306  end if
1307  if (allocated(epsm1_trcc_qbz)) then
1308    ABI_FREE(epsm1_trcc_qbz)
1309  end if
1310  if (allocated(epsm1_tmp)) then
1311    ABI_FREE(epsm1_tmp)
1312  end if
1313  if (allocated(degtab)) then
1314    ABI_FREE(degtab)
1315  end if
1316  if (allocated(ac_epsm1cqwz2)) then
1317    ABI_FREE(ac_epsm1cqwz2)
1318  end if
1319  if (allocated(ac_integr)) then
1320    ABI_FREE(ac_integr)
1321  end if
1322  if (allocated(aherm_sigc_ket)) then
1323    ABI_FREE(aherm_sigc_ket)
1324  end if
1325  if (allocated(herm_sigc_ket)) then
1326    ABI_FREE(herm_sigc_ket)
1327  end if
1328  if (Sigp%gwcomp==1) then
1329    ABI_FREE(wf1swf2_g)
1330  endif
1331  if (Sigp%gwcomp==1) then
1332    ABI_FREE(extrapolar_distrb)
1333  end if
1334  ABI_FREE(proc_distrb)
1335 
1336  call timab(431,2,tsec)
1337  call timab(424,2,tsec) ! calc_sigc_me
1338 
1339  call cwtime(cpu_time,wall_time,gflops,"stop")
1340  write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time
1341 
1342  DBG_EXIT("COLL")
1343 
1344 end subroutine calc_sigc_me

ABINIT/m_sigc [ Modules ]

[ Top ] [ Modules ]

NAME

  m_sigc

FUNCTION

  Compute matrix elements of the correlated part of the e-h self-energy

PARENTS

CHILDREN

SOURCE

15 #if defined HAVE_CONFIG_H
16 #include "config.h"
17 #endif
18 
19 #include "abi_common.h"
20 
21 module m_sigc
22 
23  use defs_basis
24  use defs_datatypes
25  use defs_abitypes
26  use m_gwdefs
27  use m_abicore
28  use m_xmpi
29  use m_xomp
30  use m_defs_ptgroups
31  use m_errors
32  use m_time
33  use m_splines
34 
35  !use m_gwdefs, only : czero_gw, cone_gw, Kron15N, Kron15W, Gau7W, &
36  !                     Kron23N, Kron23W, Gau11W, Kron31N, Kron31W, Gau15W
37  use m_hide_blas,     only : xdotc, xgemv, xgemm
38  use m_numeric_tools, only : hermitianize, imin_loc, coeffs_gausslegint
39  use m_fstrings,      only : sjoin, itoa
40  use m_geometry,      only : normv
41  use m_crystal,       only : crystal_t
42  use m_bz_mesh,       only : kmesh_t, get_BZ_item, findqg0, littlegroup_t, littlegroup_print
43  use m_gsphere,       only : gsphere_t, gsph_fft_tabs
44  use m_fft_mesh,      only : get_gftt, rotate_fft_mesh, cigfft
45  use m_vcoul,         only : vcoul_t
46  use m_wfd,           only : wfd_get_ur, wfd_t, wfd_get_cprj, wfd_barrier, wfd_change_ngfft, wfd_paw_get_aeur, &
47 &                            wfd_get_many_ur,wfd_sym_ur
48  use m_oscillators,   only : rho_tw_g, calc_wfwfg
49  use m_screening,     only : epsilonm1_results, epsm1_symmetrizer, epsm1_symmetrizer_inplace, get_epsm1
50  use m_ppmodel,       only : setup_ppmodel, ppm_get_qbz, ppmodel_t, calc_sig_ppm
51  use m_sigma,         only : sigma_t, sigma_distribute_bks
52  use m_esymm,         only : esymm_t, esymm_symmetrize_mels, esymm_failed
53  use m_pawang,        only : pawang_type
54  use m_pawtab,        only : pawtab_type
55  use m_pawfgrtab,     only : pawfgrtab_type
56  use m_pawcprj,       only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap
57  use m_pawpwij,       only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g
58  use m_paw_sym,       only : paw_symcprj
59  use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t
60 
61  implicit none
62 
63  private