TABLE OF CONTENTS
ABINIT/calc_coh_comp [ 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 coulomb singularity at q==0 (see csigme for the method used), note that in case of 3-D systems the factor 4pi in the Coulomb 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
SOURCE
1415 subroutine calc_coh_comp(iqibz,i_sz,same_band,nspinor,nsig_ab,ediff,npwc,gvec,& 1416 & ngfft,nfftot,wfg2_jk,vc_sqrt,botsq,otq,sigcohme) 1417 1418 !Arguments ------------------------------------ 1419 !scalars 1420 integer,intent(in) :: iqibz,npwc,nsig_ab,nspinor,nfftot 1421 real(dp),intent(in) :: i_sz,ediff 1422 logical,intent(in) :: same_band 1423 !arrays 1424 integer,intent(in) :: gvec(3,npwc),ngfft(18) 1425 complex(gwpc),intent(in) :: botsq(npwc,npwc),otq(npwc,npwc) 1426 complex(gwpc),intent(in) :: vc_sqrt(npwc) 1427 complex(gwpc),intent(in) :: wfg2_jk(nfftot*nsig_ab) 1428 complex(gwpc),intent(out) :: sigcohme(nsig_ab) 1429 1430 !Local variables------------------------------- 1431 !scalars 1432 integer,save :: enough=0 1433 integer :: ig,ig4,ig4x,ig4y,ig4z,igp,igmin,ispinor,ngfft1,ngfft2,ngfft3,spad,outofbox 1434 !arrays 1435 integer :: g2mg1(3) 1436 1437 ! ************************************************************************* 1438 1439 DBG_ENTER("COLL") 1440 1441 ! === Treat the case q --> 0 adequately === 1442 ! TODO Better treatment of wings 1443 igmin=1 ; if (iqibz==1) igmin=2 1444 ! 1445 ! === Partial contribution to the matrix element of Sigma_c === 1446 ! * For nspinor==2, the closure relation reads: 1447 ! $\sum_s \psi_a^*(1)\psi_b(2) = \delta_{ab} \delta(1-2)$ 1448 ! where a,b are the spinor components. As a consequence, Sigma_{COH} is always 1449 ! diagonal in spin-space and only diagonal matrix elements have to be calculated. 1450 ! MG TODO wfg2_jk should be calculated on an augmented FFT box to avoid spurious wrapping of G1-G2. 1451 ! 1452 ngfft1 = ngfft(1); ngfft2 = ngfft(2); ngfft3 = ngfft(3) 1453 sigcohme(:) = czero_gw 1454 1455 do ispinor=1,nspinor 1456 spad=(ispinor-1)*nfftot 1457 outofbox=0 1458 1459 do igp=igmin,npwc 1460 do ig=igmin,npwc 1461 1462 g2mg1 = gvec(:,igp)-gvec(:,ig) 1463 if (ANY(g2mg1(:)>ngfft(1:3)/2) .or. ANY(g2mg1(:)<-(ngfft(1:3)-1)/2)) then 1464 outofbox = outofbox+1; CYCLE 1465 end if 1466 1467 ig4x=MODULO(g2mg1(1),ngfft1) 1468 ig4y=MODULO(g2mg1(2),ngfft2) 1469 ig4z=MODULO(g2mg1(3),ngfft3) 1470 ig4= 1+ig4x+ig4y*ngfft1+ig4z*ngfft1*ngfft2 1471 1472 !MG where is neta here, ediff, otq might be close to zero depending on gwecomp 1473 sigcohme(ispinor) = sigcohme(ispinor) + & 1474 & half*wfg2_jk(spad+ig4)*vc_sqrt(ig)*vc_sqrt(igp) * botsq(ig,igp) / ( otq(ig,igp) * ( ediff -otq(ig,igp) ) ) 1475 end do 1476 end do 1477 1478 if (iqibz==1.and.same_band) then 1479 sigcohme(ispinor) = sigcohme(ispinor) + half*wfg2_jk(spad+ig4)*i_sz*botsq(1,1) / ( otq(1,1) * (ediff -otq(1,1)) ) 1480 end if 1481 end do !ispinor 1482 1483 if (outofbox/=0) then 1484 enough=enough+1 1485 if (enough<=50) then 1486 ABI_WARNING(sjoin('Number of G1-G2 pairs outside the G-sphere for Wfns: ',itoa(outofbox))) 1487 if (enough==50) then 1488 call wrtout(std_out,' ========== Stop writing Warnings ==========') 1489 end if 1490 end if 1491 end if 1492 1493 DBG_EXIT("COLL") 1494 1495 end subroutine calc_coh_comp
ABINIT/calc_sig_ppm_comp [ 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
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 !Arguments ------------------------------------ 1952 !scalars 1953 integer,intent(in) :: nomega,npwc,npwc1,npwc2,npwx,ppmodel 1954 real(dp),intent(in) :: omegame0i_io,theta_mu_minus_e0i,zcut 1955 !arrays 1956 complex(gwpc),intent(in) :: botsq(npwc,npwc1),rhotwgp(npwx),otq(npwc,npwc2) 1957 complex(gwpc),intent(inout) :: ket(npwc,nomega) 1958 1959 !Local variables------------------------------- 1960 !scalars 1961 integer :: ig,igp,io 1962 real(dp) :: den,otw,twofm1_zcut 1963 complex(gwpc) :: num,rhotwgdp_igp 1964 logical :: fully_occupied,totally_empty 1965 character(len=500) :: msg 1966 !arrays 1967 complex(gwpc),allocatable :: ket_comp(:) 1968 1969 !************************************************************************* 1970 1971 if (ppmodel/=1.and.ppmodel/=2) then 1972 write(msg,'(a,i0,a)')' The completeness trick cannot be used when ppmodel is ',ppmodel,' It should be set to 1 or 2. ' 1973 ABI_ERROR(msg) 1974 end if 1975 1976 ABI_MALLOC(ket_comp,(npwc)) 1977 ket_comp(:)=0.d0 1978 1979 fully_occupied=(abs(theta_mu_minus_e0i-1.)<0.001) 1980 totally_empty=(abs(theta_mu_minus_e0i)<0.001) 1981 1982 if(.not.(totally_empty)) then ! not totally empty 1983 twofm1_zcut=zcut 1984 do igp=1,npwc 1985 rhotwgdp_igp=rhotwgp(igp) 1986 do ig=1,npwc 1987 otw=DBLE(otq(ig,igp)) ! in principle otw -> otw - ieta 1988 num = botsq(ig,igp)*rhotwgdp_igp 1989 1990 den = omegame0i_io-otw 1991 if (den**2>zcut**2) then 1992 ket_comp(ig) = ket_comp(ig) - num/(den*otw)*theta_mu_minus_e0i 1993 end if 1994 end do !ig 1995 end do !igp 1996 end if ! not totally empty 1997 1998 if(.not.(fully_occupied)) then ! not fully occupied 1999 twofm1_zcut=-zcut 2000 2001 do igp=1,npwc 2002 rhotwgdp_igp=rhotwgp(igp) 2003 do ig=1,npwc 2004 otw=DBLE(otq(ig,igp)) ! in principle otw -> otw - ieta 2005 num = botsq(ig,igp)*rhotwgdp_igp 2006 2007 den = omegame0i_io-otw 2008 if (den**2>zcut**2) then 2009 ket_comp(ig) = ket_comp(ig) - num/(den*otw)*(1.-theta_mu_minus_e0i) 2010 end if 2011 end do !ig 2012 end do !igp 2013 end if ! not fully occupied 2014 2015 do io=1,nomega 2016 ket(:,io)=ket(:,io)+0.5*ket_comp(:) 2017 end do 2018 2019 ABI_FREE(ket_comp) 2020 2021 end subroutine calc_sig_ppm_comp
ABINIT/calc_sigc_cd [ 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.
SOURCE
1528 subroutine calc_sigc_cd(npwc,npwx,nspinor,nomega,nomegae,nomegaer,nomegaei,rhotwgp,& 1529 & omega,epsm1q,omegame0i,theta_mu_minus_e0i,ket,plasmafreq,npoles_missing,calc_poles,method) 1530 1531 !Arguments ------------------------------------ 1532 !scalars 1533 integer,intent(in) :: nomega,nomegae,nomegaei,nomegaer,npwc,npwx,nspinor 1534 integer,intent(inout) :: npoles_missing 1535 real(dp),intent(in) :: theta_mu_minus_e0i,plasmafreq 1536 !arrays 1537 real(dp),intent(in) :: omegame0i(nomega) 1538 complex(dpc),intent(in) :: omega(nomegae) 1539 complex(gwpc),intent(in) :: epsm1q(npwc,npwc,nomegae) 1540 complex(gwpc),intent(in) :: rhotwgp(npwx*nspinor) 1541 complex(gwpc),intent(inout) :: ket(nspinor*npwc,nomega) 1542 logical, intent(in), optional :: calc_poles(nomega) 1543 integer, intent(in), optional :: method 1544 1545 !Local variables------------------------------- 1546 !scalars 1547 integer, parameter :: FABIEN=1,TRAPEZOID=2,NSPLINE=3 1548 integer :: ii,ig,io,ios,ispinor,spadc,spadx,my_err,ierr,GK_LEVEL,INTMETHOD 1549 integer :: i,j 1550 real(dp) :: rt_imag,rt_real,local_one,local_zero 1551 real(dp) :: intsign,temp1,temp2,temp3,temp4 1552 real(dp) :: alph,inv_alph,beta,alphsq,betasq,inv_beta 1553 real(dp) :: re_intG,re_intK,im_intG,im_intK,GKttab,tau,ttil 1554 real(dp) :: ref,imf,r,s,r2,s2 1555 complex(dpc) :: ct,domegaleft,domegaright 1556 complex(gwpc) :: fact 1557 !arrays 1558 real(dp) :: omegame0i_tmp(nomega),tmp_x(2),tmp_y(2) 1559 real(dp) :: left(nomega),right(nomega) 1560 real(dp) :: tbeta(nomega),tinv_beta(nomega),tbetasq(nomega) 1561 real(dp) :: atermr(nomega),aterml(nomega),logup(nomega),logdown(nomega) 1562 real(dp) :: rtmp_r(nomegaer),rtmp_i(nomegaer) 1563 real(dp) :: ftab(nomegaei+2),ftab2(nomegaei+2),xtab(nomegaei+2),y(3,nomegaei+2) 1564 real(dp) :: work(nomegaei+2),work2(nomegaei+2),y2(3,nomegaei+2) 1565 complex(dpc) :: omega_imag(nomegaei+1) 1566 complex(gwpc) :: epsrho(npwc,nomegae),epsrho_imag(npwc,nomegaei+1) 1567 complex(gwpc) :: tfone(npwc,nomegaei+1),tftwo(npwc,nomegaei+1) 1568 complex(gwpc) :: weight(nomegaei+1,nomega) 1569 complex(gwpc) :: weight2(nomegaei,nomega) 1570 logical :: my_calc_poles(nomega) 1571 real(dp), allocatable :: KronN(:),KronW(:),GaussW(:),fint(:),fint2(:) 1572 !************************************************************************* 1573 1574 my_calc_poles=.TRUE.; my_err=0 1575 1576 ! Set integration method for imaginary axis 1577 INTMETHOD = FABIEN 1578 if (present(method)) then 1579 if (method==1) INTMETHOD = FABIEN 1580 if (method==2) INTMETHOD = TRAPEZOID 1581 if (method>2) then 1582 INTMETHOD = NSPLINE 1583 if (method==3) then 1584 GK_LEVEL = 15 1585 ABI_MALLOC(KronN,(GK_LEVEL)) 1586 ABI_MALLOC(KronW,(GK_LEVEL)) 1587 ABI_MALLOC(GaussW,(GK_LEVEL-8)) 1588 ABI_MALLOC(fint,(GK_LEVEL)) 1589 ABI_MALLOC(fint2,(GK_LEVEL)) 1590 KronN(:) = Kron15N(:); KronW(:) = Kron15W(:); GaussW(:) = Gau7W(:) 1591 else if (method==4) then 1592 GK_LEVEL = 23 1593 ABI_MALLOC(KronN,(GK_LEVEL)) 1594 ABI_MALLOC(KronW,(GK_LEVEL)) 1595 ABI_MALLOC(GaussW,(GK_LEVEL-12)) 1596 ABI_MALLOC(fint,(GK_LEVEL)) 1597 ABI_MALLOC(fint2,(GK_LEVEL)) 1598 KronN(:) = Kron23N(:); KronW(:) = Kron23W(:); GaussW(:) = Gau11W(:) 1599 else if (method>4) then 1600 GK_LEVEL = 31 1601 ABI_MALLOC(KronN,(GK_LEVEL)) 1602 ABI_MALLOC(KronW,(GK_LEVEL)) 1603 ABI_MALLOC(GaussW,(GK_LEVEL-16)) 1604 ABI_MALLOC(fint,(GK_LEVEL)) 1605 ABI_MALLOC(fint2,(GK_LEVEL)) 1606 KronN(:) = Kron31N(:); KronW(:) = Kron31W(:); GaussW(:) = Gau15W(:) 1607 end if 1608 end if 1609 end if 1610 1611 ! Avoid divergences in $\omega - \omega_s$. 1612 omegame0i_tmp(:)=omegame0i(:) 1613 do ios=1,nomega 1614 if (ABS(omegame0i_tmp(ios))<tol6) omegame0i_tmp(ios)=sign(tol6,omegame0i_tmp(ios)) 1615 end do 1616 1617 do ispinor=1,nspinor 1618 spadx=(ispinor-1)*npwx; spadc=(ispinor-1)*npwc 1619 1620 ! Calculate $ \sum_{Gp} (\epsilon^{-1}_{G Gp}(\omega)-\delta_{G Gp}) \rhotwgp(Gp) $ 1621 !$omp parallel do 1622 do io=1,nomegae 1623 call XGEMV('N',npwc,npwc,cone_gw,epsm1q(:,:,io),npwc,rhotwgp(1+spadx:),1,czero_gw,epsrho(:,io),1) 1624 end do 1625 1626 ! Integrand along the imaginary axis. 1627 epsrho_imag(:,1)=epsrho(:,1) 1628 epsrho_imag(:,2:nomegaei+1)=epsrho(:,nomegaer+1:nomegae) 1629 1630 ! Frequency mesh for integral along the imaginary axis. 1631 omega_imag(1)=omega(1) 1632 omega_imag(2:nomegaei+1)=omega(nomegaer+1:nomegae) 1633 1634 ! Original implementation -- saved here for reference during development 1635 ! === Perform integration along the imaginary axis === 1636 !do io=1,nomegaei+1 1637 ! if (io==1) then 1638 ! domegaleft = omega_imag(io) 1639 ! domegaright =(omega_imag(io+1)-omega_imag(io ))*half 1640 ! else if (io==nomegaei+1) then 1641 ! domegaleft =(omega_imag(io )-omega_imag(io-1))*half 1642 ! domegaright =(omega_imag(io )-omega_imag(io-1))*half 1643 ! else 1644 ! domegaleft =(omega_imag(io )-omega_imag(io-1))*half 1645 ! domegaright =(omega_imag(io+1)-omega_imag(io ))*half 1646 ! end if 1647 ! do ios=1,nomega 1648 ! omg2 = -AIMAG(omega_imag(io)+domegaright)/REAL(omegame0i_tmp(ios)) 1649 ! omg1 = -AIMAG(omega_imag(io)-domegaleft )/REAL(omegame0i_tmp(ios)) 1650 ! fact = ATAN(omg2)-ATAN(omg1) 1651 ! ket(spadc+1:spadc+npwc,ios)=ket(spadc+1:spadc+npwc,ios)+epsrho_imag(:,io)*fact 1652 ! end do 1653 !end do !io 1654 1655 !ket(spadc+1:spadc+npwc,:)=ket(spadc+1:spadc+npwc,:)/pi 1656 ! ---------------- end of original implementation ----------------------- 1657 1658 select case(INTMETHOD) 1659 case(FABIEN) 1660 ! Hopefully more effective implementation MS 04.08.2011 1661 ! Perform integration along imaginary axis using BLAS 1662 ! First calculate first and last weights 1663 weight(1,:) = ATAN(-half*AIMAG(omega_imag(2))/REAL(omegame0i_tmp(:))) 1664 domegaleft = (three*omega_imag(nomegaei+1)-omega_imag(nomegaei)) 1665 domegaright = (omega_imag(nomegaei+1)+omega_imag(nomegaei)) 1666 right(:) = -AIMAG(omega_imag(nomegaei+1)-omega_imag(nomegaei))*REAL(omegame0i_tmp(:)) 1667 left(:) = quarter*AIMAG(domegaleft)*AIMAG(domegaright) & 1668 & +REAL(omegame0i_tmp(:))*REAL(omegame0i_tmp(:)) 1669 do ios=1,nomega 1670 weight(nomegaei+1,ios) = ATAN(right(ios)/left(ios)) 1671 end do 1672 ! Calculate the rest of the weights 1673 do io=2,nomegaei 1674 domegaleft = (omega_imag(io )+omega_imag(io-1)) 1675 domegaright = (omega_imag(io+1)+omega_imag(io )) 1676 right(:) = -half*AIMAG(omega_imag(io+1)-omega_imag(io-1))*REAL(omegame0i_tmp(:)) 1677 left(:) = REAL(omegame0i_tmp(:))*REAL(omegame0i_tmp(:)) & 1678 & +quarter*AIMAG(domegaleft)*AIMAG(domegaright) 1679 do ios=1,nomega 1680 weight(io,ios) = ATAN(right(ios)/left(ios)) 1681 end do 1682 end do 1683 1684 ! Use BLAS call to perform matrix-matrix multiplication and accumulation 1685 fact = CMPLX(piinv,zero) 1686 1687 call xgemm('N','N',npwc,nomega,nomegaei+1,fact,epsrho_imag,npwc,& 1688 & weight,nomegaei+1,cone_gw,ket(spadc+1:spadc+npwc,:),npwc) 1689 1690 case(TRAPEZOID) 1691 ! Trapezoidal rule Transform omega coordinates 1692 alph = plasmafreq 1693 alphsq = alph*alph 1694 inv_alph = one/alph 1695 1696 xtab(1:nomegaei+1) = AIMAG(omega_imag(:))/(AIMAG(omega_imag(:)) + alph) 1697 xtab(nomegaei+2) = one 1698 1699 ! Efficient trapezoidal rule with BLAS calls 1700 tbeta(:) = REAL(omegame0i_tmp(:)) 1701 tbetasq(:) = tbeta(:)*tbeta(:) 1702 tinv_beta(:) = one/tbeta(:) 1703 1704 do io=1,nomegaei 1705 atermr(:) = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(io+1)-tbetasq(:)) 1706 aterml(:) = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(io )-tbetasq(:)) 1707 right(:) = ATAN((atermr(:)-aterml(:))/(one+atermr(:)*aterml(:))) 1708 logup(:) = ABS(((alphsq+tbetasq(:))*xtab(io+1)-two*tbetasq(:)) & 1709 & *xtab(io+1)+tbetasq(:)) 1710 logdown(:) = ABS(((alphsq+tbetasq(:))*xtab(io )-two*tbetasq(:)) & 1711 & *xtab(io )+tbetasq(:)) 1712 ! Trapezoid integration weights 1713 weight(io,:) = CMPLX(-(half*alph*tbeta(:)*LOG(logup(:)/logdown(:)) + tbetasq(:) & 1714 & *right(:))/(alphsq+tbetasq(:)),zero) 1715 weight2(io,:) = CMPLX(-right(:),zero) 1716 ! Linear interpolation coefficients for each section (sum over ig) 1717 tfone(:,io) = (epsrho_imag(:,io+1)-epsrho_imag(:,io)) & 1718 & /(xtab(io+1)-xtab(io)) 1719 tftwo(:,io) = epsrho_imag(:,io) - tfone(:,io)*xtab(io) 1720 end do 1721 1722 ! Calculate weights for asymptotic behaviour 1723 atermr(:) = alph*tinv_beta(:) 1724 aterml(:) = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(nomegaei+1)-tbetasq(:)) 1725 logup(:) = alphsq*xtab(nomegaei+1)*xtab(nomegaei+1) 1726 logdown(:) = ABS(((alphsq+tbetasq(:))*xtab(nomegaei+1)-two*tbetasq(:)) & 1727 & *xtab(nomegaei+1)+tbetasq(:)) 1728 right(:) = ATAN((atermr(:)-aterml(:))/(one+atermr(:)*aterml(:))) 1729 weight (nomegaei+1,:) = CMPLX(-(half*(alphsq*tinv_beta(:)*LOG(logdown(:)/logup(:)) & 1730 & - tbeta(:)*LOG(xtab(nomegaei+1)*xtab(nomegaei+1))) - alph*right(:)),zero) 1731 tfone(:,nomegaei+1) = -(zero-epsrho_imag(:,nomegaei+1)*AIMAG(omega_imag(nomegaei+1))) & 1732 /(one-xtab(nomegaei+1)) 1733 1734 ! Use BLAS call to perform matrix-matrix multiplication and accumulation 1735 fact = CMPLX(piinv,zero) 1736 1737 call xgemm('N','N',npwc,nomega,nomegaei+1,fact,tfone,npwc,& 1738 & weight ,nomegaei+1,cone_gw,ket(spadc+1:spadc+npwc,:),npwc) 1739 call xgemm('N','N',npwc,nomega,nomegaei ,fact,tftwo,npwc,& 1740 & weight2,nomegaei ,cone_gw,ket(spadc+1:spadc+npwc,:),npwc) 1741 1742 case(NSPLINE) 1743 ! Natural spline followed by Gauss-Kronrod 1744 ! Transform omega coordinates 1745 alph = plasmafreq 1746 alphsq = alph*alph 1747 inv_alph = one/alph 1748 1749 xtab(1:nomegaei+1) = AIMAG(omega_imag(:))/(AIMAG(omega_imag(:)) + alph) 1750 xtab(nomegaei+2) = one 1751 1752 ! Gauss-Kronrod integration of spline fit of f(t)/(1-t) in transformed space 1753 ! *** OPENMP SECTION *** Added by MS 1754 !!$OMP PARALLEL DO PRIVATE(ig,ftab,ftab2,s,s2,r,r2,y,y2,work,work2,beta,betasq,inv_beta, & 1755 !!$OMP intsign,io,ii,i,j,re_intG,re_intK,im_intG,im_intK,temp1,temp2,temp3,temp4, & 1756 !!$OMP ttil,tau,ref,fint,imf,fint2,GKttab) 1757 do ig=1,npwc 1758 ! Spline fit 1759 ftab (1:nomegaei+1) = REAL(epsrho_imag(ig,1:nomegaei+1))/(one-xtab(1:nomegaei+1)) 1760 ftab2(1:nomegaei+1) = AIMAG(epsrho_imag(ig,1:nomegaei+1))/(one-xtab(1:nomegaei+1)) 1761 ftab (nomegaei+2) = zero; ftab2(nomegaei+2) = zero 1762 ! Explicit calculation of spline coefficients 1763 s = zero; s2 = zero 1764 do i = 1, nomegaei+2-1 1765 r = ( ftab (i+1) - ftab (i) ) / ( xtab(i+1) - xtab(i) ) 1766 r2 = ( ftab2(i+1) - ftab2(i) ) / ( xtab(i+1) - xtab(i) ) 1767 y (2,i) = r - s; y2(2,i) = r2 - s2 1768 s = r; s2 = r2 1769 end do 1770 s = zero; s2 = zero 1771 r = zero; r2 = zero 1772 y(2,1) = zero; y2(2,1) = zero 1773 y(2,nomegaei+2) = zero; y2(2,nomegaei+2) = zero 1774 do i = 2, nomegaei+2-1 1775 y (2,i) = y (2,i) + r * y (2,i-1) 1776 y2(2,i) = y2(2,i) + r2 * y2(2,i-1) 1777 work (i) = two * ( xtab(i-1) - xtab(i+1) ) - r * s 1778 work2(i) = two * ( xtab(i-1) - xtab(i+1) ) - r2 * s2 1779 s = xtab(i+1) - xtab(i) 1780 s2 = s 1781 r = s / work (i) 1782 r2 = s2 / work2(i) 1783 end do 1784 do j = 2, nomegaei+2-1 1785 i = nomegaei+2+1-j 1786 y (2,i) = ( ( xtab(i+1) - xtab(i) ) * y (2,i+1) - y (2,i) ) / work (i) 1787 y2(2,i) = ( ( xtab(i+1) - xtab(i) ) * y2(2,i+1) - y2(2,i) ) / work2(i) 1788 end do 1789 do i = 1, nomegaei+2-1 1790 s = xtab(i+1) - xtab(i); s2 = s; 1791 r = y(2,i+1) - y(2,i); r2 = y2(2,i+1) - y2(2,i); 1792 y(3,i) = r / s; y2(3,i) = r2 / s2; 1793 y(2,i) = three * y(2,i); y2(2,i) = three * y2(2,i); 1794 y (1,i) = ( ftab (i+1) - ftab (i) ) / s - ( y (2,i) + r ) * s 1795 y2(1,i) = ( ftab2(i+1) - ftab2(i) ) / s2 - ( y2(2,i) + r2 ) * s2 1796 end do 1797 ! End of spline interpolation 1798 do ios=1,nomega 1799 beta = REAL(omegame0i_tmp(ios)) 1800 betasq = beta*beta 1801 inv_beta = one/beta 1802 intsign = sign(half*piinv,beta) 1803 beta = ABS(beta) 1804 io = 1; re_intG = zero; re_intK = zero; im_intG = zero; im_intK = zero 1805 do ii=1,GK_LEVEL 1806 do 1807 GKttab = two*alph*xtab(io+1)/(beta-(beta-alph)*xtab(io+1))-one 1808 if (GKttab > KronN(ii)) EXIT 1809 io = io + 1 1810 end do 1811 temp1 = half*(KronN(ii)+one) 1812 temp2 = temp1 - half 1813 temp3 = temp2*temp2 1814 temp4 = half/(temp3 + quarter) 1815 ttil = beta*temp1/(alph-(alph-beta)*temp1) 1816 tau = ttil - xtab(io) 1817 ref = ftab (io) + tau*(y (1,io)+tau*(y (2,io)+tau*y (3,io))) 1818 fint (ii) = -ref*(one-ttil)*temp4 1819 imf = ftab2(io) + tau*(y2(1,io)+tau*(y2(2,io)+tau*y2(3,io))) 1820 fint2(ii) = -imf*(one-ttil)*temp4 1821 re_intK = KronW(ii)*fint (ii) 1822 im_intK = KronW(ii)*fint2(ii) 1823 ket(spadc+ig,ios) = ket(spadc+ig,ios)+intsign*CMPLX(re_intK,im_intK) 1824 end do ! ii 1825 end do !ios 1826 end do !ig 1827 !!$OMP END PARALLEL DO 1828 1829 end select !intmethod 1830 1831 local_one = one 1832 local_zero = zero 1833 1834 ! ============================================ 1835 ! ==== Add contribution coming from poles ==== 1836 ! ============================================ 1837 ! First see if the contribution has been checked before the routine is entered 1838 if (present(calc_poles)) then 1839 my_calc_poles = calc_poles 1840 else ! Otherwise check locally if there is a contribution 1841 do ios=1,nomega 1842 if (omegame0i_tmp(ios)>tol12) then 1843 if ((local_one-theta_mu_minus_e0i)<tol12) my_calc_poles(ios) = .FALSE. 1844 end if 1845 if (omegame0i_tmp(ios)<-tol12) then 1846 if (theta_mu_minus_e0i<tol12) my_calc_poles(ios) = .FALSE. 1847 end if 1848 end do !ios 1849 end if 1850 1851 if (ANY(my_calc_poles(:))) then ! Make sure we only enter if necessary 1852 ! *** OPENMP SECTION *** Added by MS 1853 !!OMP !write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',xomp_get_num_threads() 1854 !$OMP PARALLEL SHARED(npwc,nomega,nomegaer,theta_mu_minus_e0i,spadc,local_one,local_zero, & 1855 !$OMP omega,epsrho,omegame0i_tmp,ket,my_calc_poles) & 1856 !$OMP PRIVATE(ig,ios,rtmp_r,rtmp_i,tmp_x,tmp_y,rt_real,rt_imag,ct,ierr) REDUCTION(+:my_err) 1857 !!OMP $ write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',xomp_get_num_threads() 1858 !$OMP DO 1859 do ig=1,npwc 1860 ! Prepare the spline interpolation by filling at once the arrays rtmp_r, rtmp_i 1861 call spline(DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),nomegaer,local_zero,local_zero,rtmp_r) 1862 call spline(DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),nomegaer,local_zero,local_zero,rtmp_i) 1863 ! call spline_complex( DBLE(omega(1:nomegaer)), epsrho(ig,1:nomegaer), nomegaer, zero, zero, rtmp ) 1864 1865 do ios=1,nomega 1866 if (.NOT.my_calc_poles(ios)) CYCLE 1867 1868 ! Interpolate real and imaginary part of epsrho at |omegame0i_tmp|. 1869 tmp_x(1) = ABS(omegame0i_tmp(ios)) 1870 call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),rtmp_r,1,tmp_x,tmp_y,ierr=ierr) 1871 if (ig==1.and.ispinor==1) my_err = my_err + ierr 1872 rt_real = tmp_y(1) 1873 1874 tmp_x(1) = ABS(omegame0i_tmp(ios)) 1875 call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),rtmp_i,1,tmp_x,tmp_y) 1876 rt_imag = tmp_y(1) 1877 !!call splint_complex(nomegaer,DBLE(omega(1:nomegaer)),epsrho(ig,1:nomegaer),rtmp,1,tmp_x,yfit) 1878 1879 ct=DCMPLX(rt_real,rt_imag) 1880 1881 if (omegame0i_tmp(ios)>tol12) then 1882 ket(spadc+ig,ios)=ket(spadc+ig,ios)+ct*(local_one-theta_mu_minus_e0i) 1883 end if 1884 if (omegame0i_tmp(ios)<-tol12) then 1885 ket(spadc+ig,ios)=ket(spadc+ig,ios)-ct*theta_mu_minus_e0i 1886 end if 1887 1888 end do !ios 1889 end do !ig 1890 !$OMP END DO 1891 !$OMP END PARALLEL 1892 end if ! ANY(my_calc_poles) 1893 end do !ispinor 1894 1895 npoles_missing = npoles_missing + my_err 1896 1897 if (INTMETHOD>2) then 1898 ABI_FREE(KronN) 1899 ABI_FREE(KronW) 1900 ABI_FREE(GaussW) 1901 ABI_FREE(fint) 1902 ABI_FREE(fint2) 1903 end if 1904 1905 end subroutine calc_sigc_cd
ABINIT/calc_sigc_me [ 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.
SOURCE
159 subroutine calc_sigc_me(sigmak_ibz,ikcalc,nomega_sigc,minbnd,maxbnd,& 160 & Dtset,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_Max,Gsph_c,Vcp,Kmesh,Qmesh,Ltg_k,& 161 & PPm,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,allQP_sym,& 162 & gwc_ngfft,rho_ngfft,rho_nfftot,rhor,use_aerhor,aepaw_rhor,sigcme_tmp) 163 164 !Arguments ------------------------------------ 165 !scalars 166 integer,intent(in) :: sigmak_ibz,ikcalc,rho_nfftot,nomega_sigc,minbnd,maxbnd 167 integer,intent(in) :: use_aerhor 168 type(crystal_t),intent(in) :: Cryst 169 type(ebands_t),target,intent(in) :: QP_BSt 170 type(kmesh_t),intent(in) :: Kmesh,Qmesh 171 type(vcoul_t),intent(in) :: Vcp 172 type(dataset_type),intent(in) :: Dtset 173 type(epsilonm1_results),intent(inout) :: Er 174 type(gsphere_t),intent(in) :: Gsph_Max,Gsph_c 175 type(littlegroup_t),intent(in) :: Ltg_k 176 type(ppmodel_t),intent(inout) :: PPm 177 type(Pseudopotential_type),intent(in) :: Psps 178 type(pawang_type),intent(in) :: pawang 179 type(sigparams_t),target,intent(in) :: Sigp 180 type(sigma_t),intent(in) :: Sr 181 type(wfdgw_t),target,intent(inout) :: Wfd,Wfdf 182 !arrays 183 integer,intent(in) :: gwc_ngfft(18),rho_ngfft(18) 184 real(dp),intent(in) :: rhor(rho_nfftot,Wfd%nspden) 185 real(dp),intent(in) :: aepaw_rhor(rho_nfftot,Wfd%nspden*use_aerhor) 186 complex(dpc),intent(out) :: sigcme_tmp(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Wfd%nsppol*Sigp%nsig_ab) 187 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat) 188 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw) 189 type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol) 190 type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw) 191 type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom*Psps%usepaw) 192 193 !Local variables ------------------------------ 194 !scalars 195 integer,parameter :: tim_fourdp2=2,ndat1=1 196 integer :: npw_k,iab,ib,ib1,ib2,ierr,ig,ii,iik,itim_q,i1,i2,npls 197 integer :: ik_bz,ik_ibz,io,iiw,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx 198 integer :: band,band1,band2,idle,rank,jik,jk_bz,jk_ibz,kb,nspinor 199 integer :: nomega_tot,nq_summed,ibsp,dimcprj_gw,npwc 200 integer :: spad,spadc1,spadc2,irow,my_nbks,ndegs,wtqm,wtqp,mod10 201 integer :: isym_kgw,isym_ki,gwc_mgfft,use_padfft,gwc_fftalga,gwc_nfftot,nfftf,mgfftf,use_padfftf 202 integer :: iwc,ifft 203 real(dp) :: cpu_time,wall_time,gflops 204 real(dp) :: e0i,fact_sp,theta_mu_minus_e0i,tol_empty,z2,en_high,gw_gsq,w_localmax,w_max 205 complex(dpc) :: ctmp,omegame0i2_ac,omegame0i_ac,ph_mkgwt,ph_mkt 206 logical :: iscompatibleFFT,q_is_gamma 207 character(len=500) :: msg,sigma_type 208 type(wave_t),pointer :: wave_sum, wave_jb 209 complex(gwpc),allocatable :: botsq(:,:),otq(:,:),eig(:,:) 210 !arrays 211 integer :: g0(3),spinor_padc(2,4),got(Wfd%nproc) 212 integer,allocatable :: proc_distrb(:,:,:),extrapolar_distrb(:,:,:,:),degtab(:,:,:) 213 integer,allocatable :: igfftcg0(:),gw_gfft(:,:),gw_gbound(:,:),irottb(:,:),ktabr(:,:) 214 integer,allocatable :: igfftfcg0(:),gboundf(:,:),ktabrf(:,:),npoles_missing(:) 215 real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),omegap(Er%nomega_i),omegap2(Er%nomega_i),q0(3),tsec(2),qbz(3) 216 real(dp) :: gl_knots(Er%nomega_i),gl_wts(Er%nomega_i) 217 real(dp) :: spinrot_kbz(4),spinrot_kgw(4) 218 real(dp),ABI_CONTIGUOUS pointer :: qp_ene(:,:,:),qp_occ(:,:,:) 219 real(dp),allocatable :: omegame0i(:), w_maxval(:) 220 complex(gwpc) :: sigcohme(Sigp%nsig_ab) 221 complex(gwpc),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:) 222 complex(gwpc),allocatable :: botsq_conjg_transp(:,:),ac_epsm1cqwz2(:,:,:) 223 complex(gwpc),allocatable :: epsm1_qbz(:,:,:),epsm1_trcc_qbz(:,:,:), epsm1_tmp(:,:) 224 complex(gwpc),allocatable :: sigc_ket(:,:),ket1(:,:),ket2(:,:) 225 complex(gwpc),allocatable :: herm_sigc_ket(:,:),aherm_sigc_ket(:,:) 226 complex(gwpc),allocatable :: rhotwg_ki(:,:) 227 complex(gwpc),allocatable :: sigcme2(:,:),sigcme_3(:),sigcme_new(:),sigctmp(:,:) 228 complex(gwpc),allocatable :: wfr_bdgw(:,:),ur_ibz(:),wf1swf2_g(:),usr_bz(:) 229 complex(gwpc),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:) 230 complex(gwpc),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:) 231 complex(gwpc),allocatable :: otq_transp(:,:) 232 complex(gwpc),ABI_CONTIGUOUS pointer :: cg_jb(:),cg_sum(:) 233 complex(dpc),allocatable :: sym_cme(:,:,:,:),sigc(:,:,:,:,:) 234 logical :: rank_mask(Wfd%nproc),can_symmetrize(Wfd%nsppol) 235 !logical :: me_calc_poles(Sr%nomega_r+Sr%nomega4sd) 236 type(sigijtab_t),pointer :: Sigcij_tab(:) 237 type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:) 238 type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:) 239 type(esymm_t),pointer :: QP_sym(:) 240 integer :: ilwrk, neigmax 241 integer :: neig(Er%nomega_i) 242 real(gwp) :: epsm1_ev(Sigp%npwc) 243 complex(gwpc),allocatable :: epsm1_sqrt_rhotw(:,:), rhotw_epsm1_rhotw(:,:,:) 244 245 !************************************************************************ 246 247 DBG_ENTER("COLL") 248 249 ! Initial check 250 ABI_CHECK(Sr%nomega_r==Sigp%nomegasr,"") 251 ABI_CHECK(Sr%nomega4sd==Sigp%nomegasrd,"") 252 ABI_CHECK(Sigp%npwc==Gsph_c%ng,"") 253 ABI_CHECK(Sigp%npwvec==Gsph_Max%ng,"") 254 255 mod10=MOD(Sigp%gwcalctyp,10) 256 257 call timab(424,1,tsec) ! calc_sigc_me 258 call timab(431,1,tsec) ! calc_sigc_me 259 call timab(432,1,tsec) ! Init 260 call cwtime(cpu_time,wall_time,gflops,"start") 261 262 ! Initialize MPI variables 263 qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ 264 265 ! Extract the symmetries of the bands for this k-point 266 QP_sym => allQP_sym(sigmak_ibz,1:Wfd%nsppol) 267 268 ! Index of the GW point in the BZ array, its image in IBZ and time-reversal === 269 jk_bz=Sigp%kptgw2bz(ikcalc) 270 call kmesh%get_BZ_item(jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt) 271 !%call kmesh%get_IBZ_item(jk_ibz,kibz,wtk) 272 273 ! TODO: the new version based of get_uug won't suppporte kptgw vectors that are not in 274 ! the IBZ since one should perform the rotation before entering the band loop 275 ! In the old version, the rotation was done in rho_tw_g 276 !ABI_CHECK(jik==1,"jik!=1") 277 !ABI_CHECK(isym_kgw==1,"isym_kgw!=1") 278 !ABI_CHECK((ABS(ph_mkgwt - cone) < tol12),"ph_mkgwt!") 279 280 spinrot_kgw=Cryst%spinrot(:,isym_kgw) 281 ib1=minbnd; ib2=maxbnd 282 283 write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,& 284 ' Calculating <nk|Sigma_c(omega)|nk> at k = ',kgw(:),ch10,& 285 ' bands n = from ',ib1,' to ',ib2,ch10 286 call wrtout(std_out, msg) 287 288 if ( Dtset%gwaclowrank > 0 ) then 289 ! Today we use the same number of eigenvectors irrespective to iw' 290 ! Tomorrow we might optimize this further 291 neigmax = MIN(Dtset%gwaclowrank,Sigp%npwc) 292 else 293 neigmax = Sigp%npwc 294 endif 295 296 ABI_MALLOC(w_maxval,(minbnd:maxbnd)) 297 w_maxval = zero 298 299 if ( ANY(gwc_ngfft(1:3) /= Wfd%ngfft(1:3)) ) call Wfd%change_ngfft(Cryst,Psps,gwc_ngfft) 300 gwc_mgfft = MAXVAL(gwc_ngfft(1:3)) 301 gwc_fftalga = gwc_ngfft(7)/100 !; gwc_fftalgc=MOD(gwc_ngfft(7),10) 302 303 if (Dtset%pawcross==1) mgfftf = MAXVAL(rho_ngfft(1:3)) 304 305 can_symmetrize = .FALSE. 306 if (Sigp%symsigma>0) then 307 can_symmetrize = .TRUE. 308 if (Sigp%gwcalctyp >= 20) then 309 do spin=1,Wfd%nsppol 310 can_symmetrize(spin) = .not. esymm_failed(QP_sym(spin)) 311 if (.not.can_symmetrize(spin)) then 312 write(msg,'(a,i0,4a)')& 313 " Symmetrization cannot be performed for spin: ",spin,ch10,& 314 " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg) 315 ABI_WARNING(msg) 316 end if 317 end do 318 end if 319 if (Wfd%nspinor == 2) ABI_WARNING("Symmetrization with nspinor=2 not implemented") 320 end if 321 322 ABI_UNUSED(Pawang%l_max) 323 324 ! Print type of calculation. 325 sigma_type = sigma_type_from_key(mod10) 326 call wrtout(std_out, sigma_type) 327 328 ! Set up logical flags for Sigma calculation 329 if (mod10==SIG_GW_AC.and.Sigp%gwcalctyp/=1) then 330 if(Sigp%gwcalctyp/=21) then 331 ABI_ERROR("not implemented") 332 else 333 write(msg,'(a34,i9)')'Constructing Sigma_c(iw) for k = ',ikcalc 334 call wrtout([std_out, ab_out], msg) 335 end if 336 end if 337 338 if (mod10==SIG_GW_AC.and.Sigp%gwcomp==1) then 339 ABI_ERROR("not implemented") 340 end if 341 342 if (mod10==SIG_GW_AC) then 343 write(msg,'(3a,i6,a,i6)')& 344 ' Using a low-rank formula for AC', ch10, & 345 ' Number of epsm1 eigenvectors retained: ',neigmax,' over ',Sigp%npwc 346 call wrtout(std_out, msg) 347 endif 348 349 350 ! Initialize some values 351 nspinor = Wfd%nspinor 352 npwc = Sigp%npwc 353 spinor_padc(:,:)=RESHAPE([0, 0, npwc, npwc, 0, npwc, npwc, 0], [2, 4]) 354 ABI_MALLOC(npoles_missing,(minbnd:maxbnd)) 355 npoles_missing=0 356 357 ! Normalization of theta_mu_minus_e0i 358 ! If nsppol==2, qp_occ $\in [0,1]$ 359 SELECT CASE (Wfd%nsppol) 360 CASE (1) 361 fact_sp=half; tol_empty=0.01 ! below this value the state is assumed empty 362 if (Wfd%nspinor==2) then 363 fact_sp=one; tol_empty=0.005 ! below this value the state is assumed empty 364 end if 365 CASE (2) 366 fact_sp=one; tol_empty=0.005 ! to be consistent and obtain similar results if a metallic 367 CASE DEFAULT ! spin unpolarized system is treated using nsppol==2 368 ABI_BUG('Wrong nsppol') 369 END SELECT 370 371 ! Allocate arrays used to accumulate the matrix elements of \Sigma_c over 372 ! k-points and bands. Note that for AC requires only the imaginary frequencies 373 !nomega_sigc=Sr%nomega_r+Sr%nomega4sd 374 !if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i 375 ! 376 ! === Define the G-G0 shifts for the FFT of the oscillators === 377 ! * Sigp%mG0 gives the MAX G0 component to account for umklapp. 378 ! * Note the size MAX(Sigp%npwx,npwc). 379 ! 380 ! === Precalculate the FFT index of $(R^{-1}(r-\tau))$ === 381 ! * S=\transpose R^{-1} and k_BZ = S k_IBZ 382 ! * irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk. 383 gwc_nfftot = PRODUCT(gwc_ngfft(1:3)) 384 ABI_MALLOC(irottb,(gwc_nfftot,Cryst%nsym)) 385 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwc_ngfft,irottb,iscompatibleFFT) 386 if (.not.iscompatibleFFT) then 387 msg = "FFT mesh is not compatible with symmetries. Results might be affected by large errors!" 388 ABI_WARNING(msg) 389 end if 390 391 ABI_MALLOC(ktabr,(gwc_nfftot,Kmesh%nbz)) 392 do ik_bz=1,Kmesh%nbz 393 isym=Kmesh%tabo(ik_bz) 394 do ifft=1,gwc_nfftot 395 ktabr(ifft,ik_bz)=irottb(ifft,isym) 396 end do 397 end do 398 ABI_FREE(irottb) 399 400 if (Psps%usepaw==1 .and. Dtset%pawcross==1) then 401 nfftf = PRODUCT(rho_ngfft(1:3)) 402 ABI_MALLOC(irottb,(nfftf,Cryst%nsym)) 403 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,rho_ngfft,irottb,iscompatibleFFT) 404 405 ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz)) 406 do ik_bz=1,Kmesh%nbz 407 isym=Kmesh%tabo(ik_bz) 408 do ifft=1,nfftf 409 ktabrf(ifft,ik_bz)=irottb(ifft,isym) 410 end do 411 end do 412 ABI_FREE(irottb) 413 end if 414 415 Sigcij_tab => Sigp%Sigcij_tab(ikcalc,1:Wfd%nsppol) 416 417 got=0 418 ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,Wfd%nsppol)) 419 call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,Wfd%nsppol,can_symmetrize,kgw,Sigp%mg0,& 420 & my_nbks,proc_distrb,got,global=.TRUE.) 421 422 write(msg,'(a,i0,a)')" Will sum ",my_nbks," (b,k,s) states in Sigma_c." 423 call wrtout(std_out, msg) 424 425 if (Sigp%gwcomp==1) then 426 en_high=MAXVAL(qp_ene(Sigp%nbnds,:,:)) + Sigp%gwencomp 427 write(msg,'(6a,e11.4,a)')ch10,& 428 ' Using the extrapolar approximation to accelerate convergence',ch10,& 429 ' with respect to the number of bands included',ch10,& 430 ' with extrapolar energy: ',en_high*Ha_eV,' [eV]' 431 call wrtout(std_out, msg) 432 ABI_MALLOC(wf1swf2_g,(gwc_nfftot*nspinor)) 433 end if 434 435 if (Sigp%gwcomp == 1) then 436 ! Setup of MPI table for extrapolar contributions. 437 ABI_MALLOC(extrapolar_distrb,(ib1:ib2,ib1:ib2,Kmesh%nbz,Wfd%nsppol)) 438 extrapolar_distrb = xmpi_undefined_rank 439 440 do spin=1,Wfd%nsppol 441 do ik_bz=1,Kmesh%nbz 442 if (ANY(proc_distrb(:,ik_bz,spin) /= xmpi_undefined_rank) ) then ! This BZ point will be calculated. 443 rank_mask = .FALSE. ! The set of node that will treat (k,s). 444 do band=1,Wfd%mband 445 rank = proc_distrb(band,ik_bz,spin) 446 if (rank /= xmpi_undefined_rank) rank_mask(rank+1)=.TRUE. 447 end do 448 do band2=ib1,ib2 449 do irow=1,Sigcij_tab(spin)%col(band2)%size1 ! Looping over the non-zero elements of sigma_ij. 450 band1 = Sigcij_tab(spin)%col(band2)%bidx(irow) 451 idle = imin_loc(got,mask=rank_mask) 452 got(idle) = got(idle)+1 453 extrapolar_distrb(band1,band2,ik_bz,spin) = idle-1 454 end do 455 end do 456 end if 457 end do 458 end do 459 460 write(msg,'(a,i0,a)')" Will treat ",COUNT(extrapolar_distrb==Wfd%my_rank)," extrapolar terms." 461 call wrtout(std_out, msg) 462 end if 463 464 ABI_MALLOC(rhotwg_ki, (npwc*nspinor, minbnd:maxbnd)) 465 rhotwg_ki=czero_gw 466 ABI_MALLOC(rhotwg, (npwc*nspinor)) 467 ABI_MALLOC(rhotwgp, (npwc*nspinor)) 468 ABI_MALLOC(vc_sqrt_qbz, (npwc)) 469 470 if (Er%mqmem == 0) then ! Use out-of-core solution for epsilon. 471 ABI_COMMENT('Reading q-slices from file. Slower but less memory.') 472 end if ! 473 474 ! Additional allocations for PAW. 475 if (Psps%usepaw==1) then 476 ABI_MALLOC(Cprj_ksum,(Cryst%natom,nspinor)) 477 call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm) 478 ! 479 ! For the extrapolar method we need the onsite terms of the PW in the FT mesh. 480 ! * gw_gfft is the set of plane waves in the FFT Box for the oscillators. 481 if (Sigp%gwcomp==1) then 482 ABI_MALLOC(gw_gfft,(3,gwc_nfftot)) 483 q0=zero 484 call get_gftt(gwc_ngfft,q0,Cryst%gmet,gw_gsq,gw_gfft) 485 ABI_MALLOC(Pwij_fft,(Psps%ntypat)) 486 call pawpwij_init(Pwij_fft,gwc_nfftot, [zero,zero,zero], gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 487 end if 488 end if ! usepaw==1 489 490 if (mod10==SIG_GW_AC) then 491 ! Calculate Gauss-Legendre quadrature knots and weights for analytic continuation 492 493 ABI_MALLOC(rhotw_epsm1_rhotw, (minbnd:maxbnd, minbnd:maxbnd, Er%nomega_i)) 494 call coeffs_gausslegint(zero,one,gl_knots,gl_wts,Er%nomega_i) 495 496 ierr = 0 497 do io=1,Er%nomega_i ! First frequencies are always real 498 if (ABS(AIMAG(one*Er%omega(Er%nomega_r+io))-(one/gl_knots(io)-one)) > 0.0001) then 499 ierr = ierr + 1 500 if (Wfd%my_rank == Wfd%master) then 501 if (io == 1) write(std_out, "(a)")"omega_file, gauss_legendre_omega (ev)" 502 write(std_out,*)io, AIMAG(Er%omega(Er%nomega_r+io)) * Ha_eV, (one/gl_knots(io)-one) * Ha_eV 503 end if 504 end if 505 end do 506 if (ierr /= 0) then 507 write(std_out, *)"Er%nomega_r:", Er%nomega_r, "Er%nomega_i:", Er%nomega_i 508 write(msg,'(3a)')& 509 'Frequencies in the SCR file are not compatible with the analytic continuation. ',ch10,& 510 'Verify the frequencies in the SCR file. ' 511 ABI_ERROR(msg) 512 end if 513 514 ! To calculate \int_0^\infty domegap f(omegap), we calculate \int_0^1 dz f(1/z-1)/z^2. 515 omegap(:)=one/gl_knots(:)-one 516 omegap2(:)=omegap(:)*omegap(:) 517 ABI_MALLOC(ac_epsm1cqwz2, (npwc, npwc, Er%nomega_i)) 518 end if 519 520 ! Calculate total number of frequencies and allocate related arrays. 521 ! sigcme2 is used to accumulate the diagonal matrix elements over k-points and 522 ! GW bands, used only in case of ppmodel 3 and 4 (TODO save memory) 523 nomega_tot=Sr%nomega_r+Sr%nomega4sd 524 ABI_MALLOC(sigcme2,(nomega_tot,ib1:ib2)) 525 ABI_MALLOC(sigcme_3,(nomega_tot)) 526 sigcme2=czero_gw; sigcme_3=czero_gw 527 528 ABI_MALLOC(sigctmp,(nomega_sigc,Sigp%nsig_ab)) 529 sigctmp=czero_gw 530 if (mod10/=SIG_GW_AC) then 531 ABI_MALLOC(sigc_ket, (npwc*nspinor, nomega_sigc)) 532 endif 533 534 #if 0 535 !TODO gmatteo: these arrays are never used in practice. Should we remove them? 536 ! Arrays storing the contribution given by the Hermitian/anti-Hermitian part of \Sigma_c 537 ABI_MALLOC(aherm_sigc_ket, (npwc*nspinor, nomega_sigc)) 538 ABI_MALLOC( herm_sigc_ket, (npwc*nspinor, nomega_sigc)) 539 #endif 540 541 sigcme_tmp=czero 542 543 ABI_MALLOC(sigc,(2,nomega_sigc,ib1:ib2,ib1:ib2,Wfd%nsppol*Sigp%nsig_ab)) 544 sigc=czero 545 546 if( mod10==SIG_QPGW_PPM .or. mod10==SIG_QPGW_CD ) then 547 ABI_MALLOC(ket1, (npwc*nspinor, nomega_tot)) 548 ABI_MALLOC(ket2, (npwc*nspinor, nomega_tot)) 549 endif 550 ABI_MALLOC(omegame0i,(nomega_tot)) 551 552 ! Here we divide the states where the QP energies are required into degenerate groups 553 ! Note however that this approach is not based on group theory, and it might lead to 554 ! spurious results in case of accidental degeneracies. 555 nq_summed=Kmesh%nbz 556 if (Sigp%symsigma > 0) then 557 call Ltg_k%print(std_out, Dtset%prtvol, mode_paral='COLL') 558 nq_summed=SUM(Ltg_k%ibzq(:)) 559 ! 560 ! Find number of degenerate subspaces and number of bands in each subspace 561 ! The tolerance is a little bit arbitrary (0.001 eV) 562 ! It could be reduced, in particular in case of nearly accidental degeneracies 563 ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,Wfd%nsppol)) 564 degtab=0 565 do spin=1,Wfd%nsppol 566 do ib=ib1,ib2 567 do jb=ib1,ib2 568 if (ABS(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin))<0.001/Ha_ev) then 569 degtab(ib,jb,spin)=1 570 end if 571 end do 572 end do 573 end do 574 end if !symsigma 575 576 write(msg,'(2a,i6,a)')ch10,' calculation status ( ',nq_summed,' to be completed):' 577 call wrtout(std_out, msg) 578 579 ! Here we have a problem in case of CD since epsm1q might be huge 580 ! TODO if single q (ex molecule) dont allocate epsm1q, avoid waste of memory 581 if ( ANY(mod10 == [SIG_GW_AC, SIG_GW_CD, SIG_QPGW_CD])) then 582 if (.not.(mod10==SIG_GW_CD.and.Er%mqmem==0)) then 583 ABI_MALLOC_OR_DIE(epsm1_qbz, (npwc, npwc, Er%nomega), ierr) 584 end if 585 end if 586 587 ! TODO epsm1_trcc_qbz is needed for SIG_GW_CD with symmetries since 588 ! the Hermitian and the anti-Hermitian part have to be symmetrized in a different way. 589 !if (mod10==SIG_QPGW_CD) then 590 if (mod10==SIG_QPGW_CD) then 591 ABI_MALLOC_OR_DIE(epsm1_trcc_qbz, (npwc, npwc, Er%nomega), ierr) 592 ABI_MALLOC(epsm1_tmp, (npwc, npwc)) 593 end if 594 595 ABI_MALLOC(igfftcg0,(Gsph_Max%ng)) 596 ABI_MALLOC(ur_ibz,(gwc_nfftot*nspinor)) 597 ABI_MALLOC(usr_bz,(gwc_nfftot*nspinor)) 598 599 if (Dtset%pawcross==1) then 600 ABI_MALLOC(igfftfcg0,(Gsph_c%ng)) 601 ABI_MALLOC(ur_ae_sum,(nfftf*nspinor)) 602 ABI_MALLOC(ur_ae_onsite_sum,(nfftf*nspinor)) 603 ABI_MALLOC(ur_ps_onsite_sum,(nfftf*nspinor)) 604 end if 605 call timab(432,2,tsec) ! Init 606 607 ! ========================================== 608 ! ==== Fat loop over k_i in the full BZ ==== 609 ! ========================================== 610 do spin=1,Wfd%nsppol 611 if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE 612 call timab(433,1,tsec) ! Init spin 613 614 ! Load wavefunctions for GW corrections 615 ! TODO: Rotate the functions here instead of calling rho_tw_g 616 ABI_MALLOC(wfr_bdgw,(gwc_nfftot*nspinor,ib1:ib2)) 617 call Wfd%get_many_ur([(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw) 618 619 if (Wfd%usepaw==1) then 620 ! Load cprj for GW states, note the indexing. 621 dimcprj_gw=nspinor*(ib2-ib1+1) 622 ABI_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1)) 623 call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm) 624 ibsp=ib1 625 do jb=ib1,ib2 626 call Wfd%get_cprj(jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 627 call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 628 call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1))) 629 ibsp=ibsp+nspinor 630 end do 631 if (Dtset%pawcross==1) then 632 ABI_MALLOC(ur_ae_bdgw,(nfftf*nspinor,ib1:ib2)) 633 ABI_MALLOC(ur_ae_onsite_bdgw,(nfftf*nspinor,ib1:ib2)) 634 ABI_MALLOC(ur_ps_onsite_bdgw,(nfftf*nspinor,ib1:ib2)) 635 do jb=ib1,ib2 636 call Wfdf%paw_get_aeur(jb,jk_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,& 637 & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum) 638 ur_ae_bdgw(:,jb)=ur_ae_sum 639 ur_ae_onsite_bdgw(:,jb)=ur_ae_onsite_sum 640 ur_ps_onsite_bdgw(:,jb)=ur_ps_onsite_sum 641 end do 642 end if 643 end if ! usepaw 644 645 call timab(433,2,tsec) ! Init spin 646 647 do ik_bz=1,Kmesh%nbz 648 ! Parallelization over k-points and spin 649 ! For the spin there is another check in the inner loop 650 if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE 651 652 call timab(434,1,tsec) ! initq 653 654 ! Find the corresponding irreducible k-point 655 call kmesh%get_BZ_item(ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt) 656 spinrot_kbz(:)=Cryst%spinrot(:,isym_ki) 657 658 ! Identify q and G0 where q + G0 = k_GW - k_i 659 kgw_m_ksum=kgw-ksum 660 call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0) 661 662 ! If symsigma, symmetrize the matrix elements. 663 ! Sum only q"s in IBZ_k. In this case elements are weighted 664 ! according to wtqp and wtqm. wtqm is for time-reversal. 665 wtqp = 1; wtqm = 0 666 if (can_symmetrize(spin)) then 667 if (Ltg_k%ibzq(iq_bz)/=1) CYCLE 668 wtqp = 0; wtqm = 0 669 do isym=1,Ltg_k%nsym_sg 670 wtqp = wtqp + Ltg_k%wtksym(1,isym,iq_bz) 671 wtqm = wtqm + Ltg_k%wtksym(2,isym,iq_bz) 672 end do 673 end if 674 675 write(msg,'(3(a,i0),a,i0)')' Sigma_c: ik_bz ',ik_bz,'/',Kmesh%nbz,", spin: ",spin,' done by mpi-rank: ',Wfd%my_rank 676 call wrtout(std_out, msg) 677 678 ! Find the corresponding irred q-point. 679 call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q) 680 q_is_gamma = normv(qbz, Cryst%gmet, "G") < GW_TOLQ0 681 682 !q_is_gamma = (normv(qbz,Cryst%gmet,"G") < 0.7) 683 !if (iq_ibz/=2.and.iq_ibz/=1) CYCLE 684 !if (ANY(qbz<=-(half-tol16)) .or. ANY(qbz>(half+tol16))) CYCLE 685 !if (q_is_gamma) then; write(std_out,*)"skipping q=Gamma"; CYCLE; end if 686 ! 687 ! Tables for the FFT of the oscillators. 688 ! a) FFT index of the G-G0. 689 ! b) gw_gbound table for the zero-padded FFT performed in rhotwg. 690 ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2)) 691 call Gsph_c%fft_tabs(g0,gwc_mgfft,gwc_ngfft,use_padfft,gw_gbound,igfftcg0) 692 693 if (ANY(gwc_fftalga == [2, 4])) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g 694 #ifdef FC_IBM 695 ! XLF does not deserve this optimization (problem with [v67mbpt][t03]) 696 use_padfft = 0 697 #endif 698 if (use_padfft==0) then 699 ABI_FREE(gw_gbound) 700 ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2*use_padfft)) 701 end if 702 703 if (Dtset%pawcross==1) then 704 ABI_MALLOC(gboundf,(2*mgfftf+8,2)) 705 call Gsph_c%fft_tabs(g0,mgfftf,rho_ngfft,use_padfftf,gboundf,igfftfcg0) 706 if ( ANY(gwc_fftalga == (/2,4/)) ) use_padfftf=0 707 if (use_padfftf==0) then 708 ABI_FREE(gboundf) 709 ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf)) 710 end if 711 end if 712 713 ! Evaluate oscillator matrix elements 714 ! $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form. 715 if (Psps%usepaw==1) then 716 ABI_MALLOC(Pwij_qg,(Psps%ntypat)) 717 q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT 718 call pawpwij_init(Pwij_qg,npwc,q0,Gsph_c%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 719 end if 720 721 if (Er%mqmem==0) then 722 ! Read q-slice of epsilon^{-1}|chi0 in Er%epsm1(:,:,:,1) (much slower but less memory). 723 call get_epsm1(Er,Vcp,0,0,Dtset%iomode,xmpi_comm_self,iqibzA=iq_ibz) 724 if (sigma_needs_ppm(Sigp)) then 725 if (Wfd%usepaw==1.and.PPm%userho==1) then 726 ! Use PAW AE rhor. 727 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,rho_nfftot,Gsph_c%gvec,& 728 rho_ngfft,aepaw_rhor(:,1),iqiA=iq_ibz) 729 else 730 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,rho_nfftot,Gsph_c%gvec,& 731 rho_ngfft,rhor(:,1),iqiA=iq_ibz) 732 end if 733 end if 734 end if 735 736 ! Symmetrize PPM parameters and epsm1 (q_IBZ --> q_BZ): 737 ! NOTE: 738 ! - We are not considering umklapp with G0/=0. In this case, 739 ! indeed the equation is different since we have to use G-G0. 740 ! A check, however, is performed in sigma. 741 ! - If gwcomp==1 and mod10 in [1,2,9], one needs both to set up botsq and epsm1_q 742 if (sigma_needs_ppm(Sigp)) then 743 ABI_MALLOC(botsq,(PPm%npwc,PPm%dm2_botsq)) 744 ABI_MALLOC(otq,(PPm%npwc,PPm%dm2_otq)) 745 ABI_MALLOC(eig,(PPm%dm_eig,PPm%dm_eig)) 746 call ppm_get_qbz(PPm,Gsph_c,Qmesh,iq_bz,botsq,otq,eig) 747 end if 748 749 if ( ANY(mod10 == [SIG_GW_AC,SIG_GW_CD,SIG_QPGW_CD])) then 750 ! Numerical integration or model GW with contour deformation or Analytic Continuation 751 ! TODO In case of AC we should symmetrize only the imaginary frequencies 752 if (mod10==SIG_GW_CD.and.Er%mqmem==0) then 753 ! Do in-place symmetrisation 754 call Epsm1_symmetrizer_inplace(iq_bz,Er%nomega, npwc,Er,Gsph_c,Qmesh,.TRUE.) 755 else 756 call Epsm1_symmetrizer(iq_bz,Er%nomega, npwc,Er,Gsph_c,Qmesh,.TRUE.,epsm1_qbz) 757 end if 758 759 if (mod10==SIG_GW_AC) then 760 761 call timab(444,1,tsec) ! ac_lrk_diag 762 ! Important to set to zero here for all procs since we're going to use a dirty reduction later 763 ! The reduction 'xmpi_sum' does not induce a significant performance loss in the tested systems 764 ac_epsm1cqwz2(:,:,:) = czero_gw 765 neig(:) = 0 766 do iiw=1,Er%nomega_i 767 ! Use the available MPI tasks to parallelize over iw' 768 if ( Dtset%gwpara == 2 .and. MODULO(iiw-1,Wfd%nproc) /= Wfd%my_rank ) CYCLE 769 770 ! Prepare the integration weights w_i 1/z_i^2 f(1/z_i-1).. 771 ! The first frequencies are always real, skip them. 772 z2=gl_knots(iiw)*gl_knots(iiw) 773 ac_epsm1cqwz2(:,:,iiw)= gl_wts(iiw)*epsm1_qbz(:,:,Er%nomega_r+iiw)/z2 774 775 ! (epsm1-1) has negative eigenvalues 776 ! after the diago, they will be sorted starting from the most negative 777 call xheev('V','L',npwc,ac_epsm1cqwz2(:,:,iiw),epsm1_ev) 778 779 ! Eliminate the spurious positive eigenvalues that may occur in harsh conditions 780 neig(iiw) = MIN(COUNT(epsm1_ev(:)<-1.0e-10_dp),neigmax) 781 782 do ilwrk=1,neig(iiw) 783 ac_epsm1cqwz2(:,ilwrk,iiw) = ac_epsm1cqwz2(:,ilwrk,iiw) * SQRT( -epsm1_ev(ilwrk) ) 784 end do 785 end do 786 if ( Dtset%gwpara == 2 ) then 787 call xmpi_sum(ac_epsm1cqwz2, Wfd%comm, ierr) 788 call xmpi_sum(neig, Wfd%comm, ierr) 789 endif 790 call timab(444,2,tsec) ! ac_lrk_diag 791 792 end if 793 794 if (mod10==SIG_QPGW_CD) then 795 ! For model GW we need transpose(conjg(epsm1_qbz)) === 796 do io=1,Er%nomega 797 epsm1_tmp(:,:) = GWPC_CONJG(epsm1_qbz(:,:,io)) 798 epsm1_trcc_qbz(:,:,io) = TRANSPOSE(epsm1_tmp) 799 end do 800 end if 801 end if !gwcalctyp 802 803 ! Get Fourier components of the Coulomb interaction in the BZ 804 ! In 3D systems, neglecting umklapp: vc(Sq,sG) = vc(q,G) = 4pi/|q+G|**2 805 ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S. 806 do ig=1,npwc 807 vc_sqrt_qbz(Gsph_c%rottb(ig,itim_q,isym_q)) = Vcp%vc_sqrt(ig,iq_ibz) 808 end do 809 810 call timab(434,2,tsec) ! initq 811 812 ! Sum over band 813 call timab(445,1,tsec) ! loop 814 do ib=1,Sigp%nbnds 815 816 ! Parallelism over spin 817 ! This processor has this k-point but what about spin? 818 if (proc_distrb(ib,ik_bz,spin)/=Wfd%my_rank) CYCLE 819 820 call Wfd%get_ur(ib,ik_ibz,spin,ur_ibz) 821 822 if (Psps%usepaw==1) then 823 ! Load cprj for point ksum, this spin or spinor and *THIS* band. 824 ! TODO MG I could avoid doing this but I have to exchange spin and bands ??? 825 ! For sure there is a better way to do this! 826 call Wfd%get_cprj(ib,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.) 827 call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum) 828 if (Dtset%pawcross==1) then 829 call Wfdf%paw_get_aeur(ib,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,& 830 ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum) 831 end if 832 end if 833 834 call timab(436,2,tsec) ! (1) 835 call timab(437,1,tsec) ! rho_tw_g 836 837 ! Get all <k-q,ib,s|e^{-i(q+G).r}|s,jb,k>, at once 838 do jb=ib1,ib2 839 840 call rho_tw_g(nspinor,npwc,gwc_nfftot,ndat1,gwc_ngfft,1,use_padfft,igfftcg0,gw_gbound,& 841 ur_ibz ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz, & 842 wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,& 843 nspinor,rhotwg_ki(:,jb)) 844 845 if (Psps%usepaw==1) then 846 ! Add on-site contribution, projectors are already in BZ !TODO Recheck this! 847 i2=jb; if (nspinor==2) i2=(2*jb-1) 848 spad=(nspinor-1) 849 call paw_rho_tw_g(cryst,Pwij_qg,npwc,nspinor,nspinor,Gsph_c%gvec, & 850 Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),rhotwg_ki(:,jb)) 851 852 if (Dtset%pawcross==1) then ! Add paw cross term 853 call paw_cross_rho_tw_g(nspinor,npwc,nfftf,rho_ngfft,1,use_padfftf,igfftfcg0,gboundf,& 854 ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum,iik,ktabrf(:,ik_bz),ph_mkt,spinrot_kbz,& 855 ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,& 856 nspinor,rhotwg_ki(:,jb)) 857 end if 858 end if 859 860 ! Multiply by the square root of the Coulomb term 861 ! In 3-D systems, the factor sqrt(4pi) is included) 862 do ii=1,nspinor 863 spad=(ii-1)*npwc 864 rhotwg_ki(spad+1:spad+npwc,jb) = rhotwg_ki(spad+1:spad+npwc,jb)*vc_sqrt_qbz(1:npwc) 865 end do 866 867 ! === Treat analytically the case q --> 0 === 868 ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma 869 ! while the Colulomb term is integrated out. 870 ! * In the scalar case we have nonzero contribution only if ib==jb 871 ! * For nspinor==2 evalute <ib,up|jb,up> and <ib,dwn|jb,dwn>, 872 ! impose orthonormalization since npwwfn might be < npwvec. 873 if (ik_bz==jk_bz) then 874 if (nspinor==1) then 875 rhotwg_ki(1,jb)=czero_gw 876 if (ib==jb) rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) 877 878 else 879 npw_k = Wfd%npwarr(ik_ibz) 880 rhotwg_ki(1, jb) = zero; rhotwg_ki(npwc+1, jb) = zero 881 if (ib==jb) then 882 ABI_CHECK(Wfd%get_wave_ptr(ib, ik_ibz, spin, wave_sum, msg) == 0, msg) 883 cg_sum => wave_sum%ug 884 ABI_CHECK(Wfd%get_wave_ptr(jb, jk_ibz, spin, wave_jb, msg) == 0, msg) 885 cg_jb => wave_jb%ug 886 ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1) 887 rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp) 888 ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1) 889 rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp) 890 ! PAW is missing 891 end if 892 end if 893 end if 894 end do !jb Got all matrix elements from minbnd up to maxbnd. 895 896 theta_mu_minus_e0i=fact_sp*qp_occ(ib,ik_ibz,spin) 897 898 ! Starting point to evaluate the derivative of Sigma and the Spectral function 899 e0i=qp_ene(ib,ik_ibz,spin) 900 901 ! Frequencies for the spectral function, e0i=qp_ene(ib,ik_ibz,spin) 902 ! FIXME the interval is not centered on eoi ! WHY? 903 if (Sr%nomega_r>0) omegame0i(1:Sr%nomega_r)=DBLE(Sr%omega_r(1:Sr%nomega_r))-e0i 904 905 call timab(437,2,tsec) ! rho_tw_g 906 call timab(443,1,tsec) ! ac_lrk_appl 907 908 if (mod10 == SIG_GW_AC) then 909 rhotw_epsm1_rhotw(:,:,:) = czero_gw 910 do iiw=1,Er%nomega_i 911 ABI_MALLOC(epsm1_sqrt_rhotw, (neig(iiw), minbnd:maxbnd)) 912 ! epsm1_sqrt_rhotw = SQRT(epsm1) * rho_tw 913 call xgemm('C','N',neig(iiw),maxbnd-minbnd+1,npwc,cone_gw,ac_epsm1cqwz2(:,:,iiw),npwc,& 914 rhotwg_ki,npwc,czero_gw,epsm1_sqrt_rhotw,neig(iiw)) 915 call xherk('L','C',maxbnd-minbnd+1,neig(iiw),one_gw,epsm1_sqrt_rhotw,neig(iiw),zero_gw,& 916 rhotw_epsm1_rhotw(:,:,iiw),maxbnd-minbnd+1) 917 918 ! Get the upper part of rhotw_epsm1_rhotw that is hermitian by construction 919 do jb=minbnd,maxbnd 920 do kb=jb+1,maxbnd 921 rhotw_epsm1_rhotw(jb,kb,iiw) = CONJG( rhotw_epsm1_rhotw(kb,jb,iiw) ) 922 end do 923 end do 924 ABI_FREE(epsm1_sqrt_rhotw) 925 end do 926 call timab(443,2,tsec) ! ac_lrk_appl 927 endif 928 929 do kb=ib1,ib2 930 call timab(438,1,tsec) ! (2) 931 932 ! Get frequencies $\omega$-\epsilon_in$ to evaluate $d\Sigma/dE$, note the spin 933 ! subtract e_KS since we have stored e_KS+ Delta \omega in Sr%omega4sd, not required for AC 934 do io=Sr%nomega_r+1,nomega_tot 935 omegame0i(io)=DBLE(Sr%omega4sd(kb,jk_ibz,io-Sr%nomega_r,spin))-e0i 936 end do 937 938 ! Get the ket \Sigma|\phi_{k,kb}> according to the method. 939 rhotwgp(:)=rhotwg_ki(:,kb) 940 941 SELECT CASE (mod10) 942 CASE (SIG_GW_PPM) 943 ! GW WITH Plasmon-Pole Model. 944 ! Note that ppmodel 3 or 4 work only in case of standard perturbative approach! 945 ! Moreover, for ppmodel 3 and 4, spinorial case is not allowed 946 sigc_ket = czero_gw 947 call calc_sig_ppm(PPm,nspinor,npwc,nomega_tot,rhotwgp,botsq,otq,& 948 omegame0i,Sigp%zcut,theta_mu_minus_e0i,eig,npwc,sigc_ket,sigcme_3) 949 950 if (PPm%model==3.or.PPm%model==4) then 951 sigcme2(:,kb)=sigcme2(:,kb) + (wtqp+wtqm)*DBLE(sigcme_3(:)) + (wtqp-wtqm)*j_gw*AIMAG(sigcme_3(:)) 952 end if 953 954 CASE (SIG_GW_AC) 955 ! GW with Analytic continuation. 956 ! This part is so optimized for AC that there is nothing to do here ! 957 958 CASE (SIG_GW_CD) 959 ! GW with contour deformation. 960 ! Check if pole contributions need to be summed 961 ! (this avoids unnecessary splint calls and saves time) 962 !me_calc_poles = .TRUE. 963 sigc_ket = czero_gw 964 do io=1,nomega_tot 965 if (omegame0i(io)>=zero.AND.(ABS(one-theta_mu_minus_e0i)>zero)) then 966 !me_calc_poles(io) = .TRUE. 967 if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io)) 968 else if (omegame0i(io)<zero.AND.(ABS(theta_mu_minus_e0i)>zero)) then 969 !me_calc_poles(io) = .TRUE. 970 if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io)) 971 end if 972 end do 973 974 ! Check memory saving 975 if (Er%mqmem == 0) then 976 call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 977 Er%omega,Er%epsm1(:,:,:,1),omegame0i,theta_mu_minus_e0i,sigc_ket,Dtset%ppmfrq,npoles_missing(kb),& 978 method=Dtset%cd_frqim_method) 979 else 980 call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 981 Er%omega,epsm1_qbz,omegame0i,theta_mu_minus_e0i,sigc_ket,Dtset%ppmfrq,npoles_missing(kb),& 982 method=Dtset%cd_frqim_method) 983 end if 984 985 #if 0 986 if (wtqm/=0) then 987 call calc_sigc_cd(npwc,npwc,nspinor,,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 988 Er%omega,epsm1_trcc_qbz,omegame0i,theta_mu_minus_e0i,aherm_sigc_ket,Dtset%ppmfrq,npoles_missing(kb),& 989 method=Dtset%cd_frqim_method) 990 991 herm_sigc_ket = half*(sigc_ket + aherm_sigc_ket) 992 aherm_sigc_ket = half*(sigc_ket - aherm_sigc_ket) 993 else 994 herm_sigc_ket = sigc_ket 995 aherm_sigc_ket = czero_gw 996 end if 997 #endif 998 999 CASE (SIG_QPGW_PPM) 1000 ! MODEL GW calculation WITH PPm TODO Spinor not tested. 1001 ! Calculate \Sigma(E_k) |k> to obtain <j|\Sigma(E_k)|k> 1002 ABI_MALLOC(sigcme_new,(nomega_tot)) 1003 sigc_ket = czero_gw 1004 ket1 = czero_gw 1005 ket2 = czero_gw 1006 1007 call calc_sig_ppm(PPm,nspinor,npwc,nomega_tot,rhotwgp,botsq,otq,& 1008 omegame0i,Sigp%zcut,theta_mu_minus_e0i,eig,npwc,ket1,sigcme_new) 1009 1010 if (Sigp%gwcalctyp==28) then 1011 if (PPm%model/=1.and.PPm%model/=2) then 1012 ! This is needed to have npwc=PPm%dm2_botsq=PPm%dm2_otq 1013 write(msg,'(3a)')& 1014 'For the time being, gwcalctyp=28 cannot be used with ppmodel=3,4.',ch10,& 1015 'Use another Plasmon Pole Model when gwcalctyp=28.' 1016 ABI_ERROR(msg) 1017 end if 1018 ABI_MALLOC(botsq_conjg_transp,(PPm%dm2_botsq,npwc)) 1019 botsq_conjg_transp=TRANSPOSE(botsq) ! Keep these two lines separated, otherwise gfortran messes up 1020 botsq_conjg_transp=CONJG(botsq_conjg_transp) 1021 ABI_MALLOC(otq_transp,(PPm%dm2_otq,PPm%npwc)) 1022 otq_transp=TRANSPOSE(otq) 1023 1024 call calc_sig_ppm(PPm,nspinor,npwc,nomega_tot,rhotwgp,botsq_conjg_transp,otq_transp,& 1025 omegame0i,Sigp%zcut,theta_mu_minus_e0i,eig,npwc,ket2,sigcme_3) 1026 1027 ABI_FREE(botsq_conjg_transp) 1028 ABI_FREE(otq_transp) 1029 sigc_ket= half*(ket1+ket2) 1030 else 1031 sigc_ket= ket1 1032 end if 1033 1034 ABI_FREE(sigcme_new) 1035 1036 CASE (SIG_QPGW_CD) 1037 ! MODEL GW with numerical integration. 1038 ! Check if pole contributions need to be summed 1039 ! (this avoids unnecessary splint calls and saves time) 1040 !me_calc_poles = .TRUE. 1041 sigc_ket = czero_gw 1042 ket1 = czero_gw 1043 ket2 = czero_gw 1044 1045 do io=1,nomega_tot 1046 if (omegame0i(io)>=zero.AND.(ABS(one-theta_mu_minus_e0i)>zero)) then 1047 !me_calc_poles(io) = .TRUE. 1048 if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io)) 1049 else if (omegame0i(io)<zero.AND.(ABS(theta_mu_minus_e0i)>zero)) then 1050 !me_calc_poles(io) = .TRUE. 1051 if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io)) 1052 end if 1053 end do 1054 1055 ! Calculate \Sigma(E_k)|k> to obtain <j|\Sigma(E_k)|k> 1056 call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 1057 Er%omega,epsm1_qbz,omegame0i,theta_mu_minus_e0i,ket1,Dtset%ppmfrq,npoles_missing(kb),& 1058 method=Dtset%cd_frqim_method) 1059 1060 if (Sigp%gwcalctyp==29) then 1061 ! Calculate \Sigma^*(E_k)|k> to obtain <k|\Sigma(E_k)|j>^* 1062 call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,Er%nomega,Er%nomega_r,Er%nomega_i,rhotwgp,& 1063 Er%omega,epsm1_trcc_qbz,omegame0i,theta_mu_minus_e0i,ket2,Dtset%ppmfrq,npoles_missing(kb),& 1064 method=Dtset%cd_frqim_method) 1065 sigc_ket = half*(ket1+ket2) 1066 else 1067 sigc_ket = ket1 1068 end if 1069 1070 CASE DEFAULT 1071 ABI_ERROR(sjoin("Unsupported value for mod10:", itoa(mod10))) 1072 END SELECT 1073 1074 if (Sigp%gwcomp==1) then 1075 ! TODO spinor not implemented 1076 call calc_sig_ppm_comp(npwc,nomega_tot,rhotwgp,botsq,otq,DBLE(Sr%egw(kb,jk_ibz,spin)-en_high),& 1077 Sigp%zcut,theta_mu_minus_e0i,sigc_ket,PPm%model,npwc,PPm%dm2_botsq,PPm%dm2_otq) 1078 end if 1079 1080 call timab(438,2,tsec) ! 1081 call timab(439,1,tsec) ! sigma_me 1082 1083 ! Loop over the non-zero row elements of this column. 1084 ! 1) If gwcalctyp<20 : only diagonal elements since QP==KS. 1085 ! 2) If gwcalctyp>=20: only off-diagonal elements connecting states with same character. 1086 do irow=1,Sigcij_tab(spin)%col(kb)%size1 1087 jb = Sigcij_tab(spin)%col(kb)%bidx(irow) 1088 rhotwg=rhotwg_ki(:,jb) 1089 1090 ! Calculate <\phi_j|\Sigma_c|\phi_k> 1091 ! Different freqs according to method (AC or Perturbative), see nomega_sigc. 1092 if (mod10==SIG_GW_AC) then 1093 sigctmp(:,:) = czero_gw 1094 do iab=1,Sigp%nsig_ab 1095 do io=1,nomega_sigc 1096 omegame0i_ac = Sr%omega_i(io)-qp_ene(ib,ik_ibz,spin) 1097 omegame0i2_ac = omegame0i_ac*omegame0i_ac 1098 do iiw=1,Er%nomega_i 1099 sigctmp(io,iab) = sigctmp(io,iab) + piinv * rhotw_epsm1_rhotw(jb,kb,iiw) * omegame0i_ac & 1100 / (omegame0i2_ac + omegap2(iiw)) 1101 end do 1102 end do 1103 end do 1104 else 1105 do iab=1,Sigp%nsig_ab 1106 spadc1 = spinor_padc(1, iab); spadc2 = spinor_padc(2, iab) 1107 do io=1,nomega_sigc 1108 sigctmp(io,iab) = XDOTC(npwc,rhotwg(spadc1+1:),1,sigc_ket(spadc2+1:,io),1) 1109 end do 1110 end do 1111 end if 1112 1113 if (Sigp%gwcomp==1) then 1114 ! Evaluate Extrapolar term TODO this does not work with spinor 1115 if (extrapolar_distrb(jb,kb,ik_bz,spin) == Wfd%my_rank ) then 1116 ! Do it once as it does not depend on the ib index summed over. 1117 extrapolar_distrb(jb,kb,ik_bz,spin) = xmpi_undefined_rank 1118 #if 1 1119 call calc_wfwfg(ktabr(:,jk_ibz),jik, spinrot_kgw, & ! TODO: why jk_ibz? 1120 gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g) 1121 #else 1122 call calc_wfwfg(ktabr(:,jk_bz),jik,spinrot_kgw,& 1123 gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g) 1124 #endif 1125 1126 if (Psps%usepaw==1) then 1127 i1=jb; i2=kb 1128 if (nspinor==2) then 1129 i1=(2*jb-1); i2=(2*kb-1) 1130 end if 1131 spad=(nspinor-1) 1132 call paw_rho_tw_g(cryst, Pwij_fft,gwc_nfftot,Sigp%nsig_ab,nspinor, & 1133 gw_gfft,Cprj_kgw(:,i1:i1+spad),Cprj_kgw(:,i2:i2+spad),wf1swf2_g) 1134 1135 if (Dtset%pawcross==1) then ! Add paw cross term 1136 call paw_cross_rho_tw_g(nspinor,npwc,nfftf,rho_ngfft,1,use_padfftf,igfftfcg0,gboundf,& 1137 ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,& 1138 ur_ae_bdgw(:,kb),ur_ae_onsite_bdgw(:,kb),ur_ps_onsite_bdgw(:,kb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,& 1139 nspinor,wf1swf2_g) 1140 end if 1141 end if 1142 1143 ! The static contribution from completeness relation is calculated once === 1144 call calc_coh_comp(iq_ibz,Vcp%i_sz,(jb==kb),nspinor,Sigp%nsig_ab,DBLE(Sr%egw(kb,jk_ibz,spin)-en_high),& 1145 npwc,Gsph_c%gvec,gwc_ngfft,gwc_nfftot,wf1swf2_g,vc_sqrt_qbz,botsq,otq,sigcohme) 1146 1147 do io=1,nomega_sigc 1148 sigctmp(io,:) = sigctmp(io,:)+sigcohme(:) 1149 end do 1150 end if ! gwcomp==1 1151 end if ! gwcom==1 1152 1153 ! Accumulate and, in case, symmetrize matrix elements of Sigma_c 1154 do iab=1,Sigp%nsig_ab 1155 is_idx=spin; if (nspinor==2) is_idx=iab 1156 1157 sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + & 1158 (wtqp+wtqm)*DBLE(sigctmp(:,iab)) + (wtqp-wtqm)*j_gw*AIMAG(sigctmp(:,iab)) 1159 1160 sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + wtqp* sigctmp(:,iab) 1161 sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + wtqm*CONJG(sigctmp(:,iab)) 1162 ! TODO this should be the contribution coming from the anti-hermitian part. 1163 end do 1164 end do ! irow used to calculate matrix elements of $\Sigma$ 1165 1166 ! shaltaf (030406): this has to be done in a clean way later. 1167 ! TODO does not work with spinor. 1168 if (mod10==SIG_GW_PPM.and.(PPm%model==3.or.PPm%model==4)) then 1169 sigcme_tmp(:,kb,kb,spin)= sigcme2(:,kb) 1170 sigc(1,:,kb,kb,spin)= sigcme2(:,kb) 1171 sigc(2,:,kb,kb,spin)= czero 1172 end if 1173 1174 call timab(439,2,tsec) ! csigme(SigC) 1175 end do !kb to calculate matrix elements of $\Sigma$ 1176 end do !ib 1177 1178 call timab(445,2,tsec) ! csigme(SigC) 1179 1180 ! Deallocate k-dependent quantities. 1181 ABI_FREE(gw_gbound) 1182 if (Dtset%pawcross==1) then 1183 ABI_FREE(gboundf) 1184 end if 1185 1186 if (sigma_needs_ppm(Sigp)) then 1187 ABI_FREE(botsq) 1188 ABI_FREE(otq) 1189 ABI_FREE(eig) 1190 end if 1191 if (Psps%usepaw==1) then 1192 call pawpwij_free(Pwij_qg) 1193 ABI_FREE(Pwij_qg) 1194 end if 1195 1196 end do ! ik_bz 1197 1198 ABI_FREE(wfr_bdgw) 1199 if (Wfd%usepaw==1) then 1200 call pawcprj_free(Cprj_kgw) 1201 ABI_FREE(Cprj_kgw) 1202 if (Dtset%pawcross==1) then 1203 ABI_FREE(ur_ae_bdgw) 1204 ABI_FREE(ur_ae_onsite_bdgw) 1205 ABI_FREE(ur_ps_onsite_bdgw) 1206 end if 1207 end if 1208 end do ! spin 1209 1210 ABI_FREE(sigcme2) 1211 ABI_FREE(sigcme_3) 1212 ABI_FREE(igfftcg0) 1213 if (Dtset%pawcross==1) then 1214 ABI_FREE(igfftfcg0) 1215 end if 1216 1217 ! Gather contributions from all the CPUs 1218 call timab(440,1,tsec) ! wfd_barrier 1219 call timab(440,2,tsec) ! wfd_barrier 1220 call timab(441,1,tsec) ! xmpi_sum 1221 1222 call xmpi_sum(sigcme_tmp, Wfd%comm, ierr) 1223 call xmpi_sum(sigc, Wfd%comm, ierr) 1224 call timab(441,2,tsec) ! xmpi_sum 1225 1226 ! Multiply by constants. In 3D systems sqrt(4pi) is included in vc_sqrt_qbz === 1227 sigcme_tmp = sigcme_tmp /(Cryst%ucvol*Kmesh%nbz) 1228 sigc = sigc /(Cryst%ucvol*Kmesh%nbz) 1229 1230 ! If we have summed over the IBZ_q now we have to average over degenerate states 1231 ! Presently only diagonal terms are considered 1232 ! TODO QP-SCGW required a more involved approach, there is a check in sigma 1233 ! TODO it does not work if nspinor==2. 1234 call timab(442,1,tsec) ! final ops 1235 1236 do spin=1,Wfd%nsppol 1237 if (can_symmetrize(spin)) then 1238 if (mod10==SIG_GW_AC) then ! FIXME here there is a problem in case of AC with symmetries 1239 ABI_MALLOC(sym_cme, (Sr%nomega_i, ib1:ib2, ib1:ib2, Sigp%nsig_ab)) 1240 else 1241 ABI_MALLOC(sym_cme, (nomega_tot, ib1:ib2, ib1:ib2, Sigp%nsig_ab)) 1242 end if 1243 sym_cme=czero 1244 1245 ! Average over degenerate diagonal elements 1246 ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results. 1247 ! another good reason to use a strict criterion for the tolerance on eigenvalues. 1248 do ib=ib1,ib2 1249 ndegs=0 1250 do jb=ib1,ib2 1251 if (degtab(ib,jb,spin)==1) then 1252 if (nspinor == 1) then 1253 sym_cme(:, ib, ib, 1) = sym_cme(:, ib, ib, 1) + SUM(sigc(:,:,jb,jb,spin), DIM=1) 1254 else 1255 do ii=1,Sigp%nsig_ab 1256 sym_cme(:, ib, ib, ii) = sym_cme(:, ib, ib, ii) + SUM(sigc(:,:,jb,jb,ii), dim=1) 1257 end do 1258 end if 1259 end if 1260 ndegs = ndegs + degtab(ib,jb,spin) 1261 end do 1262 sym_cme(:,ib,ib,:) = sym_cme(:,ib,ib,:) / ndegs 1263 end do 1264 1265 if (Sigp%gwcalctyp >= 20) then 1266 do iwc=1,nomega_sigc 1267 call esymm_symmetrize_mels(QP_sym(spin),ib1,ib2,sigc(:,iwc,:,:,spin),sym_cme(iwc,:,:,1)) 1268 end do 1269 end if 1270 1271 ! Copy symmetrized values 1272 do ib=ib1,ib2 1273 do jb=ib1,ib2 1274 !if (mod10==SIG_GW_AC.and.average_real) CYCLE ! this is to check another scheme in case of AC 1275 if (nspinor == 1) then 1276 sigcme_tmp(:,ib,jb,spin) = sym_cme(:,ib,jb,1) 1277 else 1278 do ii=1,Sigp%nsig_ab 1279 sigcme_tmp(:,ib,jb,ii) = sym_cme(:,ib,jb,ii) 1280 end do 1281 end if 1282 end do 1283 end do 1284 ABI_FREE(sym_cme) 1285 end if 1286 end do 1287 1288 ! Reconstruct the full sigma matrix from the upper triangle (only for HF, SEX and COHSEX) 1289 !if (Sigp%gwcalctyp>=20 .and. sigma_is_herm(Sigp) ) then 1290 ! ABI_CHECK(nspinor==1,"cannot hermitianize non-collinear sigma!") 1291 ! do spin=1,Wfd%nsppol 1292 ! do io=1,nomega_sigc 1293 ! call hermitianize(sigcme_tmp(io,:,:,spin),"Upper") 1294 ! end do 1295 ! end do 1296 !end if 1297 1298 ! GW with contour deformation: check on the number of poles not included. 1299 if (ANY(mod10 == [SIG_GW_CD,SIG_QPGW_CD])) then 1300 call xmpi_sum(npoles_missing, Wfd%comm, ierr) 1301 npls = SUM(npoles_missing) 1302 if (npls>0) then 1303 ABI_WARNING(sjoin("Total number of missing poles for contour deformation method:", itoa(npls))) 1304 do band=minbnd,maxbnd 1305 npls = npoles_missing(band) 1306 if (npls>0) then 1307 write(msg,'(a,2(i0,a))')" For band ",band," there are ",npls," missing poles" 1308 call wrtout(std_out, msg) 1309 end if 1310 end do 1311 end if 1312 ! Print data on the maximum value needed for the screening along the real axis 1313 w_localmax = MAXVAL(w_maxval) 1314 call xmpi_max(w_localmax,w_max, Wfd%comm, ierr) 1315 write(msg,'(a,f12.5,a)') ' Max omega value used in W(omega): ',w_max*Ha_eV,' [eV]' 1316 call wrtout(std_out, msg) 1317 end if 1318 call timab(442,2,tsec) ! final ops 1319 1320 ! =========================== 1321 ! ==== Deallocate memory ==== 1322 ! =========================== 1323 if (Psps%usepaw==1) then 1324 ABI_SFREE(gw_gfft) 1325 call pawcprj_free(Cprj_ksum) 1326 ABI_FREE(Cprj_ksum) 1327 if (allocated(Pwij_fft)) then 1328 call pawpwij_free(Pwij_fft) 1329 ABI_FREE(Pwij_fft) 1330 end if 1331 if (Dtset%pawcross==1) then 1332 ABI_FREE(ur_ae_sum) 1333 ABI_FREE(ur_ae_onsite_sum) 1334 ABI_FREE(ur_ps_onsite_sum) 1335 ABI_FREE(ktabrf) 1336 end if 1337 end if 1338 1339 ABI_FREE(npoles_missing) 1340 ABI_FREE(ur_ibz) 1341 ABI_FREE(usr_bz) 1342 ABI_FREE(ktabr) 1343 ABI_FREE(rhotwg_ki) 1344 ABI_FREE(rhotwg) 1345 ABI_FREE(rhotwgp) 1346 ABI_FREE(vc_sqrt_qbz) 1347 ABI_FREE(omegame0i) 1348 ABI_FREE(sigctmp) 1349 ABI_FREE(sigc) 1350 ABI_FREE(w_maxval) 1351 1352 ABI_SFREE(sigc_ket) 1353 ABI_SFREE(ket1) 1354 ABI_SFREE(ket2) 1355 ABI_SFREE(epsm1_qbz) 1356 ABI_SFREE(epsm1_trcc_qbz) 1357 ABI_SFREE(epsm1_tmp) 1358 ABI_SFREE(degtab) 1359 ABI_SFREE(ac_epsm1cqwz2) 1360 ABI_SFREE(rhotw_epsm1_rhotw) 1361 ABI_SFREE(aherm_sigc_ket) 1362 ABI_SFREE(herm_sigc_ket) 1363 ABI_SFREE(wf1swf2_g) 1364 ABI_SFREE(extrapolar_distrb) 1365 ABI_SFREE(proc_distrb) 1366 1367 call timab(431,2,tsec) 1368 call timab(424,2,tsec) ! calc_sigc_me 1369 1370 call cwtime(cpu_time,wall_time,gflops,"stop") 1371 write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time 1372 1373 DBG_EXIT("COLL") 1374 1375 end subroutine calc_sigc_me
ABINIT/m_sigc [ Modules ]
NAME
m_sigc
FUNCTION
Compute matrix elements of the correlated part of the e-h self-energy
SOURCE
10 #if defined HAVE_CONFIG_H 11 #include "config.h" 12 #endif 13 14 #include "abi_common.h" 15 16 module m_sigc 17 18 use defs_basis 19 use m_gwdefs 20 use m_abicore 21 use m_xmpi 22 use m_xomp 23 use m_defs_ptgroups 24 use m_errors 25 use m_time 26 use m_splines 27 use m_dtset 28 29 use defs_datatypes, only : pseudopotential_type, ebands_t 30 use m_hide_blas, only : xdotc, xgemv, xgemm, xherk 31 use m_hide_lapack, only : xheev 32 use m_numeric_tools, only : hermitianize, imin_loc, coeffs_gausslegint 33 use m_fstrings, only : sjoin, itoa 34 use m_geometry, only : normv 35 use m_crystal, only : crystal_t 36 use m_bz_mesh, only : kmesh_t, findqg0, littlegroup_t 37 use m_gsphere, only : gsphere_t 38 use m_fft_mesh, only : get_gftt, rotate_fft_mesh, cigfft 39 use m_vcoul, only : vcoul_t 40 use m_wfd, only : wfdgw_t, wave_t 41 use m_oscillators, only : rho_tw_g, calc_wfwfg 42 use m_screening, only : epsilonm1_results, epsm1_symmetrizer, epsm1_symmetrizer_inplace, get_epsm1 43 use m_ppmodel, only : setup_ppmodel, ppm_get_qbz, ppmodel_t, calc_sig_ppm 44 use m_sigma, only : sigma_t, sigma_distribute_bks 45 use m_esymm, only : esymm_t, esymm_symmetrize_mels, esymm_failed 46 use m_pawang, only : pawang_type 47 use m_pawtab, only : pawtab_type 48 use m_pawfgrtab, only : pawfgrtab_type 49 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap 50 use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g 51 use m_paw_sym, only : paw_symcprj 52 use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t 53 use m_hide_lapack, only : xheev 54 55 implicit none 56 57 private