TABLE OF CONTENTS
- ABINIT/m_chi0tk
- m_chi0tk/accumulate_chi0_q0
- m_chi0tk/accumulate_chi0sumrule
- m_chi0tk/accumulate_head_wings_imagw
- m_chi0tk/accumulate_sfchi0_q0
- m_chi0tk/approxdelta
- m_chi0tk/assemblychi0_sym
- m_chi0tk/assemblychi0sf
- m_chi0tk/calc_kkweight
- m_chi0tk/chi0_bbp_mask
- m_chi0tk/completechi0_deltapart
- m_chi0tk/hilbert_transform
- m_chi0tk/hilbert_transform_headwings
- m_chi0tk/make_transitions
- m_chi0tk/mkrhotwg_sigma
- m_chi0tk/output_chi0sumrule
- m_chi0tk/setup_spectral
- m_chi0tk/symmetrize_afm_chi0tk
ABINIT/m_chi0tk [ Modules ]
NAME
m_chi0tk
FUNCTION
This module provides tools for the computation of the irreducible polarizability.
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (MG, FB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_chi0tk 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 use m_xomp 29 use m_sort 30 31 use defs_datatypes, only : ebands_t 32 use m_gwdefs, only : GW_TOL_DOCC, czero_gw, cone_gw, one_gw, em1params_t, j_gw 33 use m_fstrings, only : sjoin, itoa 34 use m_hide_blas,only : xgerc, xgemm, xherk, xher 35 use m_crystal, only : crystal_t 36 use m_gsphere, only : gsphere_t 37 use m_bz_mesh, only : littlegroup_t, kmesh_t 38 use m_wfd, only : wfdgw_t 39 40 implicit none 41 42 private 43 44 public :: assemblychi0_sym 45 public :: symmetrize_afm_chi0 46 public :: accumulate_chi0_q0 47 public :: accumulate_head_wings_imagw 48 public :: accumulate_sfchi0_q0 49 public :: assemblychi0sf 50 public :: approxdelta 51 public :: setup_spectral 52 public :: hilbert_transform 53 public :: hilbert_transform_headwings 54 public :: completechi0_deltapart 55 public :: output_chi0sumrule 56 public :: accumulate_chi0sumrule 57 public :: make_transitions 58 public :: chi0_bbp_mask
m_chi0tk/accumulate_chi0_q0 [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
accumulate_chi0_q0
FUNCTION
Update the independent particle susceptibility at q==0 for the contribution of one pair of occupied-unoccupied band, for each frequency. This routine takes advantage of the symmetries of the little group of the external q-point to symmetrize the contribution arising from the input k-point located in the IBZ_q. It computes: $ \chi_0(G1,G2,io) = \chi_0(G1,G2,io)+\sum_S (rhotwg(G1)*rhotwg^\dagger(G2))*green_w(io) $ where S is a symmetry in reciprocal space. The matrix elements of the gradient operator and [V_{nl},r] are symmetrized as well.
INPUTS
ik_bz=Index of the k-point whose contribution has to be added to chi0. isym_kbz=Index of the symmetry such that k_bz = IS k_ibz itim_kbz=2 if time-reversal has to be used to obtain k_bz, 1 otherwise. npwepG0=Maximum number of G vectors rhotwg(npwepG0)=Oscillator matrix elements corresponding to an occupied-unoccupied pair of states. rhotwx(3,nspinor**2)=Matrix element of the operator $-i[H,r]/(e1-e2) = -i r$ in reciprocal lattice units. green_w(nomega)=Frequency dependent part of the Green function. Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point. %timrev=2 it time-reversal is used, 1 otherwise. %nsym_sg=Number of space group symmetries. %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point. Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors %ng=number of G vectors in the enlarged sphere, actually MUST be equal to the size of rhotwg. %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion. %phmGt(ng,nsym)=phase factors associated to non-simmorphic operations. Cryst<crystal_t>=Structure defining the unit cell and its symmetries %nsym=Number of symmetries. %symrec(3,3,nsym)=Symmetry operation in reciprocal space (reduced coordinates) Ep<em1params_t>=Parameters of the chi0 calculation. %npwe=number of plane waves in chi0. %symchi=1 if symmetrization has to be performed. %nomega=number of frequencies in chi0.
SIDE EFFECTS
chi0(npwe,npwe,nomega)= Updated independent-particle susceptibility matrix in reciprocal space at q==0. chi0_head(3,3,Ep%nomega)=Head. chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)=Lower wing. chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)=Upper wing.
NOTES
1) Symmetrization of the oscilator matrix elements.
If Sq = q then M_G( Sk,q)= e^{-i(q+G).\tau} M_{ S^-1G} (k,q)
If -Sq = q then M_G(-Sk,q)= e^{-i(q+G).\tau} M_{-S^-1G}^*(k,q)
In the case of umklapps:
If Sq = q+G0 then M_G( Sk,q)= e^{-i(q+G).\tau} M_{ S^-1(G-G0} (k,q)
If -Sq = q+G0 then M_G(-Sk,q)= e^{-i(q+G).\tau} M_{-S^-1(G-G0)}^*(k,q)
In the equation below there is no need to take into account the phases due to q.t
as they cancel each other in the scalar product ==> only phmGt(G,isym)=e^{-iG.\tau} is needed.
2) Symmetrization of the matrix elements of the position operator.
<Sk,b|\vec r| Sk,b'> = <k b| R\vec r + \tau|k b'>
where S is one of the symrec operation, R and \tau is the corresponding
operation in real space. The term involving the fractional translation is zero provided that b /= b'.
SOURCE
617 subroutine accumulate_chi0_q0(is_metallic,ik_bz,isym_kbz,itim_kbz,gwcomp,nspinor,npwepG0,Ep,Cryst,Ltg_q,Gsph_epsG0,& 618 chi0,rhotwx,rhotwg,green_w,green_enhigh_w,deltaf_b1b2,chi0_head,chi0_lwing,chi0_uwing) 619 620 !Arguments ------------------------------------ 621 !scalars 622 logical,intent(in) :: is_metallic 623 integer,intent(in) :: ik_bz,isym_kbz,itim_kbz,npwepG0,nspinor,gwcomp 624 real(dp),intent(in) :: deltaf_b1b2 625 type(littlegroup_t),intent(in) :: Ltg_q 626 type(gsphere_t),target,intent(in) :: Gsph_epsG0 627 type(crystal_t),intent(in) :: Cryst 628 type(em1params_t),intent(in) :: Ep 629 !arrays 630 complex(gwpc),intent(in) :: rhotwg(npwepG0) 631 complex(gwpc),intent(in) :: rhotwx(3, nspinor**2) 632 complex(gwpc),intent(inout) :: chi0(Ep%npwe*Ep%nI, Ep%npwe*Ep%nJ, Ep%nomega) 633 complex(dpc),intent(in) :: green_w(Ep%nomega), green_enhigh_w(Ep%nomega) 634 complex(dpc),intent(inout) :: chi0_head(3, 3, Ep%nomega) 635 complex(dpc),intent(inout) :: chi0_lwing(Ep%npwe*Ep%nI, Ep%nomega, 3) 636 complex(dpc),intent(inout) :: chi0_uwing(Ep%npwe*Ep%nJ, Ep%nomega, 3) 637 638 !Local variables------------------------------- 639 !scalars 640 integer :: itim,io,isym,idir,jdir,isymop,nsymop,npwe,nomega 641 real(gwp) :: dr 642 complex(gwpc) :: dd 643 !character(len=500) :: msg 644 !arrays 645 integer,ABI_CONTIGUOUS pointer :: Sm1G(:) 646 complex(dpc) :: mir_kbz(3) 647 complex(gwpc),allocatable :: rhotwg_sym(:,:) 648 complex(gwpc), ABI_CONTIGUOUS pointer :: phmGt(:) 649 650 !************************************************************************ 651 652 ABI_UNUSED(deltaf_b1b2) 653 654 npwe = ep%npwe; nomega = ep%nomega 655 656 select case (Ep%symchi) 657 case (0) 658 ! Do not use symmetries. 659 ! Symmetrize rhotwg in the full BZ and accumulate over the full BZ i.e. 660 ! chi0(G1,G2,io) = chi0(G1,G2,io) + (rhotwg(G1)*conjg(rhotwg(G2)))*green_w(io) 661 ! 662 ! The non-analytic term is symmetrized for this k-point in the BZ according to: 663 ! rhotwg(1) = S^-1q * rhotwx_ibz 664 ! rhotwg(1) = -S^-1q * conjg(rhotwx_ibz) if time-reversal is used. 665 666 ! Multiply elements G1,G2 of rhotwg by green_w(io) and accumulate in chi0(G1,G2,io) 667 668 !$OMP PARALLEL DO PRIVATE(dr, dd) IF (nomega > 2) 669 do io=1,nomega 670 ! Check if green_w(io) is real (=> pure imaginary omega) 671 ! and that it is not a metal 672 ! then the corresponding chi0(io) is hermitian 673 if (ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic) then 674 dr = green_w(io) 675 call xher('U', npwe, dr, rhotwg, 1, chi0(:,:,io), npwe) 676 else 677 dd = green_w(io) 678 call xgerc(npwe, npwe, dd, rhotwg, 1, rhotwg, 1, chi0(:,:,io), npwe) 679 endif 680 end do 681 682 ! === Accumulate heads and wings for each small q === 683 ! FIXME extrapolar method should be checked!! 684 ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'} 685 if (nspinor == 1) then 686 mir_kbz = (3-2*itim_kbz) * matmul(Cryst%symrec(:,:,isym_kbz), rhotwx(:,1)) 687 else 688 mir_kbz = (3-2*itim_kbz) * matmul(Cryst%symrec(:,:,isym_kbz), sum(rhotwx(:,1:2), dim=2)) 689 end if 690 if (itim_kbz == 2) mir_kbz = conjg(mir_kbz) 691 692 ! here we might take advantage of Hermiticity along Im axis in RPA (see mkG0w) 693 do idir=1,3 694 do io=1,nomega 695 chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_w(io) * mir_kbz(idir)*conjg(rhotwg) 696 chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_w(io) * rhotwg*conjg(mir_kbz(idir)) 697 if (gwcomp == 1) then 698 ! Add contribution due to extrapolar technique. 699 chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_enhigh_w(io) * mir_kbz(idir)*conjg(rhotwg) 700 chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_enhigh_w(io) * rhotwg*conjg(mir_kbz(idir)) 701 end if 702 end do 703 end do 704 705 ! Accumulate the head. 706 do io=1,nomega 707 do jdir=1,3 708 do idir=1,3 709 chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_w(io) * mir_kbz(idir)*conjg(mir_kbz(jdir)) 710 if (gwcomp == 1) then 711 ! Add contribution due to extrapolar technique. 712 chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_enhigh_w(io) * mir_kbz(idir)*conjg(mir_kbz(jdir)) 713 end if 714 end do 715 end do 716 end do 717 718 case (1) 719 ! Use symmetries to reconstruct the integrand. 720 721 nsymop = count(Ltg_q%wtksym(:,:,ik_bz) == 1) 722 ABI_MALLOC(rhotwg_sym, (npwe, nsymop)) 723 isymop = 0 724 725 ! Loop over symmetries of the space group and time-reversal. 726 do isym=1,Ltg_q%nsym_sg 727 do itim=1,Ltg_q%timrev 728 729 if (Ltg_q%wtksym(itim, isym, ik_bz) == 1) then 730 ! This operation belongs to the little group and has to be considered to reconstruct the BZ. 731 phmGt => Gsph_epsG0%phmGt (1:npwe, isym) ! In the 2 lines below note the slicing (1:npwe) 732 Sm1G => Gsph_epsG0%rottbm1(1:npwe, itim, isym) 733 734 isymop = isymop + 1 735 select case (itim) 736 case (1) 737 rhotwg_sym(1:npwe, isymop) = rhotwg(Sm1G(1:npwe)) * phmGt(1:npwe) 738 case (2) 739 rhotwg_sym(1:npwe,isymop) = conjg(rhotwg(Sm1G(1:npwe))) * phmGt(1:npwe) 740 case default 741 ABI_BUG(sjoin('Wrong value of itim:', itoa(itim))) 742 end select 743 744 ! === Accumulate heads and wings for each small q === 745 ! FIXME extrapolar method should be checked!! 746 747 ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'} 748 if (nspinor == 1) then 749 mir_kbz = (3-2*itim) * matmul(Cryst%symrec(:,:,isym), rhotwx(:,1)) 750 else 751 mir_kbz = (3-2*itim) * matmul(Cryst%symrec(:,:,isym), sum(rhotwx(:,1:2), dim=2)) 752 end if 753 754 if (itim == 2) mir_kbz = conjg(mir_kbz) 755 756 ! here we might take advantage of Hermiticity along Im axis in RPA (see mkG0w) 757 do idir=1,3 758 do io=1,nomega 759 chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_w(io) * mir_kbz(idir) * conjg(rhotwg_sym(:,isymop)) 760 chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_w(io) * rhotwg_sym(:,isymop) * conjg(mir_kbz(idir)) 761 if (gwcomp == 1) then 762 ! Add contribution due to extrapolar technique. 763 chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_enhigh_w(io) * mir_kbz(idir) * conjg(rhotwg_sym(:,isymop)) 764 chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_enhigh_w(io) * rhotwg_sym(:,isymop) * conjg(mir_kbz(idir)) 765 end if 766 end do 767 end do 768 769 ! Accumulate the head. 770 do io=1,nomega 771 do jdir=1,3 772 do idir=1,3 773 chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_w(io) * mir_kbz(idir) * conjg(mir_kbz(jdir)) 774 if (gwcomp == 1) then 775 ! Add contribution due to extrapolar technique. 776 chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_enhigh_w(io)*mir_kbz(idir) * conjg(mir_kbz(jdir)) 777 end if 778 end do 779 end do 780 end do 781 782 end if !wtksym 783 end do !itim 784 end do !isym 785 786 ! Multiply rhotwg_sym by green_w(io) and accumulate in chi0(G,Gp,io) 787 !$OMP PARALLEL DO PRIVATE(dr, dd) IF (nomega > 2) 788 do io=1,nomega 789 ! Check if green_w(io) is real (=> pure imaginary omega) and that it is not a metal 790 ! then the corresponding chi0(io) is hermitian 791 if (ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic) then 792 dr = green_w(io) 793 call xherk('U', 'N', npwe, nsymop, dr, rhotwg_sym, npwe, one_gw, chi0(:,:,io), npwe) 794 else 795 dd = green_w(io) 796 call xgemm('N', 'C', npwe, npwe, nsymop, dd, rhotwg_sym, npwe, rhotwg_sym, npwe, cone_gw, chi0(:,:,io), npwe) 797 endif 798 end do 799 800 ABI_FREE(rhotwg_sym) 801 802 case default 803 ABI_BUG(sjoin('Wrong value of symchi:', itoa(Ep%symchi))) 804 end select 805 806 end subroutine accumulate_chi0_q0
m_chi0tk/accumulate_chi0sumrule [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
accumulate_chi0sumrule
FUNCTION
Accumulate the contribution to the sum rule for Im chi0 arising from a single transition. Eventually symmetrize it using the symmetry operations of the little group of q.
INPUTS
ik_bz=Index of the k-point in the full BZ whose contribution has to be added and symmetrized. symchi=1 if we are summing over IBZ_q and symmetrization has to be performed. npwe=Number of G vectors in chi0. npwepG0=Number of G vectors in the "enlarged" sphere to treat umklapp. factor=factor entering the expression. delta_ene=Transition energy. Ltg_q=<littlegroup_t>= Structure gathering information on the little group of the external q. Gsph_epsG0=<gsphere_t>=Info on the G-sphere for chi0. rhotwg(npwepG0)=Fouriet transform of u_{b1 k-q} u_{b2 k} in the "enlarged" sphere.
OUTPUT
See SIDES EFFECTS SIDES EFFECTS chi0sumrule(npwe)= In input the sum rule calculated so far, In output the contribution of this transition is accounted for, and, eventually, symmetrized. using the symmetry operations of the little group of the external q.
SOURCE
2103 subroutine accumulate_chi0sumrule(ik_bz,symchi,npwe,factor,delta_ene,& 2104 & Ltg_q,Gsph_epsG0,npwepG0,rhotwg,chi0sumrule) 2105 2106 !Arguments ------------------------------------ 2107 !scalars 2108 integer,intent(in) :: ik_bz,npwe,npwepG0,symchi 2109 real(dp),intent(in) :: delta_ene,factor 2110 type(gsphere_t),intent(in) :: Gsph_epsG0 2111 type(littlegroup_t),target,intent(in) :: Ltg_q 2112 !arrays 2113 real(dp),intent(inout) :: chi0sumrule(npwe) 2114 complex(gwpc),intent(in) :: rhotwg(npwepG0) 2115 2116 !Local variables------------------------------- 2117 !scalars 2118 integer :: isym,itim 2119 !character(len=500) :: msg 2120 !arrays 2121 integer,allocatable :: Sm1_gmG0(:) 2122 integer, ABI_CONTIGUOUS pointer :: gmG0(:) 2123 complex(gwpc),allocatable :: rhotwg_sym(:) 2124 2125 !************************************************************************ 2126 2127 ! Accumulating the sum rule on chi0. 2128 ! Eq.(5.284) in G. D. Mahan Many-Particle Physics 3rd edition [[cite:Mahan2000]] 2129 2130 SELECT CASE (symchi) 2131 CASE (0) 2132 ! Do not use symmetries, sum is performed in the full BZ. 2133 chi0sumrule(:)=chi0sumrule(:) + factor*delta_ene*ABS(rhotwg(1:npwe))**2 2134 2135 CASE (1) 2136 ! Symmetrize the contribution in the full BZ. 2137 ABI_MALLOC(rhotwg_sym,(npwe)) 2138 ABI_MALLOC(Sm1_gmG0,(npwe)) 2139 2140 do itim=1,Ltg_q%timrev 2141 do isym=1,Ltg_q%nsym_sg 2142 if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then 2143 ! This operation belongs to the little group and has to be used to reconstruct the BZ === 2144 ! In the following 2 lines mind the slicing (1:npwe) 2145 gmG0 => Ltg_q%igmG0(1:npwe,itim,isym) 2146 Sm1_gmG0(1:npwe)=Gsph_epsG0%rottbm1(gmG0(1:npwe),itim,isym) 2147 rhotwg_sym(1:npwe)=rhotwg(Sm1_gmG0) 2148 2149 chi0sumrule(:)=chi0sumrule(:) + factor*delta_ene*ABS(rhotwg_sym(1:npwe))**2 2150 end if 2151 end do !isym 2152 end do !itim 2153 2154 ABI_FREE(rhotwg_sym) 2155 ABI_FREE(Sm1_gmG0) 2156 2157 CASE DEFAULT 2158 ABI_BUG(sjoin('Wrong value for symchi:', itoa(symchi))) 2159 END SELECT 2160 2161 end subroutine accumulate_chi0sumrule
m_chi0tk/accumulate_head_wings_imagw [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
accumulate_head_wings_imagw
FUNCTION
m_chi0tk/accumulate_sfchi0_q0 [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
accumulate_sfchi0_q0
FUNCTION
Update the spectral function of the independent particle susceptibility at q==0 for the contribution of one pair of occupied-unoccupied band, for each frequency. If symchi==1, the symmetries belonging to the little group of the external point q are used to reconstruct the contributions in the full Brillouin zone. In this case, the equation implented is: $ chi0(G1,G2,io)=chi0(G1,G2,io)+\sum_S (rhotwg(G1)*rhotwg^\dagger(G2))* \delta(\omega -trans) $ where S is a symmetry belonging to the little group of q. The subroutine also performs the symmetrization of the matrix elements of the gradient operator and of the commutator [V_{nl},r] with the position operator.
INPUTS
ikbz=Index in the BZ of the k-point whose contribution to chi0 has to be added, if we use symmetries, the contribution to chi0 by this k-point has to be symmetrized. isym_kbz=Index of the symmetry such as k_bz = IS k_ibz itim_kbz=2 if time-reversal has to be used to obtain k_bz, 1 otherwise. my_wl,my_wr=min and Max frequency index treated by this processor. npwe=Number of plane waves used to describe chi0. npwepG0=Maximum number of G vectors to account for umklapps. nomega=Number of frequencies in the imaginary part. rhotwg(npwepG0)=Oscillator matrix elements corresponding to an occupied-unoccupied pair of states. rhotwx(3,nspinor**2)=Matrix elements of the gradient and of the commutator of the non-local operator with the position operator. The second term is present only if inclvkb=1,2. Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors %ng=number of G vectors in the enlarged sphere. It MUST be equal to the size of rhotwg %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion %phmGt(ng,nsym)=phase factors associated to non-symmorphic operations Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point. %timrev=2 it time-reversal is used, 1 otherwise %nsym_sg=Number of space group symmetries %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point %flag_umklp(timrev,nsym)= flag for umklapp processes if 1 that the particular operation (IS) requires a G_o to preserve Q, 0 otherwise Cryst<crystal_t>=Info on unit cell and it symmetries %nsym=Number of symmetry operations. %symrec(3,3,nsym)=Symmetry operations in reciprocal space (reduced coordinates).
OUTPUT
(see side effects)
SIDE EFFECTS
sf_chi0(npwe,npwe,my_wl:my_wr)=Updated spectral function at q==0. sf_lwing(npwe,my_wl:my_wr,3)=Updated lower wing of the spectral function. sf_uwing(npwe,mw_wl:my_wr,3)=Updated upper wing of the spectral function. sf_head(3,3,my_wl:my_wr)=Updated head of the spectral function.
SOURCE
1013 subroutine accumulate_sfchi0_q0(ikbz,isym_kbz,itim_kbz,nspinor,symchi,npwepG0,npwe,Cryst,Ltg_q,Gsph_epsG0,& 1014 & factocc,my_wl,iomegal,wl,my_wr,iomegar,wr,rhotwx,rhotwg,nomegasf,sf_chi0,sf_head,sf_lwing,sf_uwing) 1015 1016 !Arguments ------------------------------------ 1017 !scalars 1018 integer,intent(in) :: ikbz,my_wl,my_wr,nomegasf,npwe,npwepG0,nspinor 1019 integer,intent(in) :: isym_kbz,itim_kbz,symchi,iomegal,iomegar 1020 real(dp),intent(in) :: factocc,wl,wr 1021 type(littlegroup_t),intent(in) :: Ltg_q 1022 type(gsphere_t),target,intent(in) :: Gsph_epsG0 1023 type(crystal_t),intent(in) :: Cryst 1024 !arrays 1025 complex(gwpc),intent(in) :: rhotwg(npwepG0) 1026 complex(gwpc),intent(in) :: rhotwx(3,nspinor**2) 1027 complex(gwpc),intent(inout) :: sf_chi0(npwe,npwe,my_wl:my_wr) 1028 complex(dpc),intent(inout) :: sf_head(3,3,my_wl:my_wr) 1029 complex(dpc),intent(inout) :: sf_lwing(npwe,my_wl:my_wr,3) 1030 complex(dpc),intent(inout) :: sf_uwing(npwe,my_wl:my_wr,3) 1031 1032 !Local variables------------------------------- 1033 !scalars 1034 integer :: itim,isym,idir,jdir 1035 complex(gwpc) :: num 1036 character(len=500) :: msg 1037 !arrays 1038 integer, ABI_CONTIGUOUS pointer :: Sm1G(:) 1039 complex(dpc) :: mir_kbz(3) 1040 complex(gwpc), ABI_CONTIGUOUS pointer :: phmGt(:) 1041 complex(gwpc),allocatable :: rhotwg_sym(:) 1042 1043 !************************************************************************ 1044 1045 if (iomegal<my_wl .or. iomegar>my_wr) then 1046 write(msg,'(3a,2(a,i0,a,i0,a))')ch10,& 1047 & 'Indices out of boundary ',ch10,& 1048 & ' my_wl = ',my_wl,' iomegal = ',iomegal,ch10,& 1049 & ' my_wr = ',my_wr,' iomegar = ',iomegar,ch10 1050 ABI_BUG(msg) 1051 end if 1052 1053 SELECT CASE (symchi) 1054 CASE (0) 1055 ! 1056 ! Calculation without symmetries 1057 ! rhotwg(1)= R^-1q*rhotwx_ibz 1058 ! rhotwg(1)=-R^-1q*conjg(rhotwx_ibz) for inversion 1059 if (wl<huge(zero)*1.d-11) then 1060 !this is awful but it is still a first coding 1061 ! Num is single precision needed for cgerc check factocc 1062 num=-wl*factocc 1063 call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,sf_chi0(:,:,iomegal),npwe) 1064 end if 1065 ! Last point, must accumulate left point but not the right one 1066 if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then 1067 num=-wr*factocc 1068 call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,sf_chi0(:,:,iomegar),npwe) 1069 end if 1070 1071 ! ================================ 1072 ! ==== Update heads and wings ==== 1073 ! ================================ 1074 1075 ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'} 1076 if (nspinor == 1) then 1077 mir_kbz =(3-2*itim_kbz) * MATMUL(Cryst%symrec(:,:,isym_kbz),rhotwx(:,1)) 1078 else 1079 mir_kbz = (3-2*itim_kbz) * MATMUL(Cryst%symrec(:,:,isym_kbz), sum(rhotwx(:,1:2), dim=2)) 1080 end if 1081 if (itim_kbz==2) mir_kbz=conjg(mir_kbz) 1082 1083 do jdir=1,3 1084 if (wl<huge(zero)*1.d-11) then 1085 ! this is awful but it is still a first coding 1086 ! Num is single precision needed for cgerc check factocc 1087 num=-wl*factocc 1088 sf_uwing(:,iomegal,jdir) = sf_uwing(:,iomegal,jdir) + num * mir_kbz(jdir) * conjg(rhotwg(1:npwepG0)) 1089 sf_lwing(:,iomegal,jdir) = sf_lwing(:,iomegal,jdir) + num * rhotwg(1:npwepG0) * conjg(mir_kbz(jdir)) 1090 do idir=1,3 1091 sf_head(idir,jdir,iomegal) = sf_head(idir,jdir,iomegal) + num * mir_kbz(idir) * conjg(mir_kbz(jdir)) 1092 end do 1093 end if 1094 1095 ! Last point, must accumulate left point but not the right one 1096 if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then 1097 num=-wr*factocc 1098 sf_uwing(:,iomegar,jdir) = sf_uwing(:,iomegar,jdir) + num * mir_kbz(jdir) * conjg(rhotwg(1:npwepG0)) 1099 sf_lwing(:,iomegar,jdir) = sf_lwing(:,iomegar,jdir) + num * rhotwg(1:npwepG0) * conjg(mir_kbz(jdir)) 1100 do idir=1,3 1101 sf_head(idir,jdir,iomegar) = sf_head(idir,jdir,iomegar) + num * mir_kbz(idir) * conjg(mir_kbz(jdir)) 1102 end do 1103 end if 1104 end do ! jdir 1105 1106 CASE (1) 1107 ! === Notes on the symmetrization of oscillator matrix elements === 1108 ! If Sq = q then M_G( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1G} (k,q) 1109 ! If -Sq = q then M_G(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1G}^*(k,q) 1110 ! 1111 ! In case of an umklapp process 1112 ! If Sq = q+G_o then M_G( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1(G-G_o} (k,q) 1113 ! If -Sq = q+G_o then M_G(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1(G-G-o)}^*(k,q) 1114 ! 1115 ! rhotwg(1)= R^-1q*rhotwx_ibz 1116 ! rhotwg(1)=-R^-1q*conjg(rhotwx_ibz) for inversion 1117 1118 ABI_MALLOC(rhotwg_sym,(npwe)) 1119 1120 ! Loop over symmetries of the space group and time-reversal 1121 do isym=1,Ltg_q%nsym_sg 1122 do itim=1,Ltg_q%timrev 1123 1124 if (Ltg_q%wtksym(itim,isym,ikbz)==1) then 1125 ! This operation belongs to the little group and has to be considered to reconstruct the BZ 1126 phmGt => Gsph_epsG0%phmGt(1:npwe,isym) ! In these 2 lines mind the slicing (1:npwe) 1127 Sm1G => Gsph_epsG0%rottbm1(1:npwe,itim,isym) 1128 1129 SELECT CASE (itim) 1130 CASE (1) 1131 rhotwg_sym(1:npwe)=rhotwg(Sm1G(1:npwe))*phmGt(1:npwe) 1132 CASE (2) 1133 rhotwg_sym(1:npwe)=conjg(rhotwg(Sm1G(1:npwe)))*phmGt(1:npwe) 1134 CASE DEFAULT 1135 ABI_BUG(sjoin('Wrong value of itim:', itoa(itim))) 1136 END SELECT 1137 1138 ! Multiply elements G,Gp of rhotwg_sym*num and accumulate in sf_chi0(G,Gp,io) 1139 if (wl<huge(zero)*1.d-11) then 1140 num=-wl*factocc 1141 call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,sf_chi0(:,:,iomegal),npwe) 1142 end if 1143 1144 ! Last point, must accumulate left point but not the right one 1145 if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then 1146 num=-wr*factocc 1147 call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,sf_chi0(:,:,iomegar),npwe) 1148 end if 1149 1150 ! Accumulate heads and wings. 1151 ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'} 1152 if (nspinor == 1) then 1153 mir_kbz =(3-2*itim) * MATMUL(Cryst%symrec(:,:,isym),rhotwx(:,1)) 1154 else 1155 mir_kbz = (3-2*itim) * MATMUL(Cryst%symrec(:,:,isym), sum(rhotwx(:,1:2), dim=2)) 1156 end if 1157 if (itim==2) mir_kbz=conjg(mir_kbz) 1158 1159 do jdir=1,3 1160 if (wl<huge(zero)*1.d-11) then 1161 ! this is awful but it is still a first coding 1162 ! Num is single precision needed for cgerc check factocc 1163 num=-wl*factocc 1164 sf_uwing(:,iomegal,jdir) = sf_uwing(:,iomegal,jdir) + num * mir_kbz(jdir) * conjg(rhotwg_sym(1:npwe)) 1165 sf_lwing(:,iomegal,jdir) = sf_lwing(:,iomegal,jdir) + num * rhotwg_sym(1:npwe) * conjg(mir_kbz(jdir)) 1166 do idir=1,3 1167 sf_head(idir,jdir,iomegal) = sf_head(idir,jdir,iomegal) + num * mir_kbz(idir) * conjg(mir_kbz(jdir)) 1168 end do 1169 end if 1170 1171 ! Last point, must accumulate left point but not the right one 1172 if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then 1173 num=-wr*factocc 1174 sf_uwing(:,iomegar,jdir) = sf_uwing(:,iomegar,jdir) + num * mir_kbz(jdir) * conjg(rhotwg_sym(1:npwe)) 1175 sf_lwing(:,iomegar,jdir) = sf_lwing(:,iomegar,jdir) + num * rhotwg_sym(1:npwe) * conjg(mir_kbz(jdir)) 1176 do idir=1,3 1177 sf_head(idir,jdir,iomegar) = sf_head(idir,jdir,iomegar) + num * mir_kbz(idir) * conjg(mir_kbz(jdir)) 1178 end do 1179 end if 1180 end do ! jdir 1181 1182 end if !wtksym 1183 end do !inv 1184 end do !isym 1185 ABI_FREE(rhotwg_sym) 1186 1187 CASE DEFAULT 1188 ABI_BUG(sjoin('Wrong value of symchi:', itoa(symchi))) 1189 END SELECT 1190 1191 end subroutine accumulate_sfchi0_q0
m_chi0tk/approxdelta [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
approxdelta
FUNCTION
Approximate the Dirac function using two methods: method 1) a triangular funtion centered at the value egwdiff_re, Eq 17 of PRB 74, 035101 (2006) [[cite:Shishkin2006]] method 2) a gaussian of witdth ep%spsmear expandended in Taylor series (at the moment only the 0-th moments) Subroutine needed to implement the calculation of the polarizability using the spectral representation as proposed in: PRB 74, 035101 (2006) [[cite:Shishkin2006]] and PRB 61, 7172 (2000) [[cite:Miyake2000]]
INPUTS
nomegasf=number of frequencies in the grid for Im \chi_0 omegasf(0:nomega+1)= frequencies (real) egwdiff_re = transition energy where the delta function is centered method= 1: a triangular shaped function used to approximated the delta 2: gaussian approximation with standard deviation (smear) smear= used only in case of method==2, defines the width of the gaussian
OUTPUT
wl = weight associated to omegal (last omega wich is smaller than egwdiff_re wr = weight associate to omegar (first omega larger than egwdff_re iomegal= index in the array omegasf of the last frequency < egwdiff iomegar= index in the array omegasf of the first frequency > egwdiff
SOURCE
1433 subroutine approxdelta(nomegasf,omegasf,egwdiff_re,smear,iomegal,iomegar,wl,wr,spmeth) 1434 1435 !Arguments ------------------------------------ 1436 !scalars 1437 integer,intent(in) :: nomegasf,spmeth 1438 integer,intent(out) :: iomegal,iomegar 1439 real(dp),intent(in) :: egwdiff_re,smear 1440 real(dp),intent(out) :: wl,wr 1441 !arrays 1442 real(dp),intent(in) :: omegasf(nomegasf) 1443 1444 !Local variables------------------------------- 1445 integer :: io,iomega 1446 real(dp) :: omegal,omegar,deltal,deltar 1447 !character(len=500) :: msg 1448 ! ************************************************************************* 1449 1450 iomega=-999 1451 do io=nomegasf,1,-1 1452 if (omegasf(io)<egwdiff_re) then 1453 iomega=io; EXIT 1454 end if 1455 end do 1456 1457 iomegal=iomega ; omegal=omegasf(iomegal) 1458 iomegar=iomegal+1; omegar=omegasf(iomegar) 1459 1460 SELECT CASE (spmeth) 1461 CASE (1) 1462 ! Weights for triangular shaped function 1463 wr= (egwdiff_re-omegal)/(omegar-omegal) 1464 wl= -(egwdiff_re-omegar)/(omegar-omegal) 1465 1466 CASE (2) 1467 ! Weights for gaussian method (0-th moment) 1468 deltal=(egwdiff_re-omegal)/smear 1469 deltar=(omegar-egwdiff_re)/smear 1470 if (deltar>=deltal) then 1471 wl=EXP(-deltal*deltal) 1472 ! this value is used to avoid double counting and speed-up 1473 wr=huge(one)*1.d-10 1474 else 1475 wl=huge(one)*1.d-10 1476 wr=exp(-deltal*deltal) 1477 end if 1478 1479 CASE DEFAULT 1480 ABI_BUG(sjoin('Wrong value for spmeth:', itoa(spmeth))) 1481 END SELECT 1482 1483 end subroutine approxdelta
m_chi0tk/assemblychi0_sym [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
assemblychi0_sym
FUNCTION
Update the independent particle susceptibility for the contribution of one pair of occupied-unoccupied band, for each frequency. If symchi=1 the expression is symmetrized taking into account the symmetries of the little group associated to the external q-point. Compute chi0(G1,G2,io)=chi0(G1,G2,io)+\sum_S \hat S (rhotwg(G1)*rhotwg*(G2))*green_w(io) where S are the symmetries of the little group associated to the external q-point.
INPUTS
nspinor=Number of spinorial components. ik_bz=Index of the k-point in the BZ array whose contribution has to be symmetrized and added to cchi0 npwepG0=Maximum number of G vectors taking into account possible umklapp G0, ie enlarged sphere G-G0 rhotwg(npwe)=Oscillator matrix elements for this k-point and the transition that has to be summed green_w(nomega)=frequency dependent part coming from the green function Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors %ng=number of G vectors in the enlarged sphere, actually MUST be equal to the size of rhotwg %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion %phmGt(ng,nsym)=phase factors associated to non-simmorphic operations Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point. %timrev=2 it time-reversal is used, 1 otherwise %nsym_sg=Number of space group symmetries %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point %flag_umklp(timrev,nsym)= flag for umklapp processes if 1 that the particular operation (IS) requires a G_o to preserve Q, 0 otherwise %igmG0(npwepG0,timrev,nsym) index of G-G0 in the array gvec Ep<em1params_t>=Parameters related to the calculation of chi0/epsilon^-1 %symchi %nomega=number of frequencies %npwe=number of plane waves for epsilon (input variable)
OUTPUT
(see side effects)
SIDE EFFECTS
chi0(npwe,npwe,nomega)=independent-particle susceptibility matrix in reciprocal space
SOURCE
109 subroutine assemblychi0_sym(is_metallic,ik_bz,nspinor,Ep,Ltg_q,green_w,npwepG0,rhotwg,Gsph_epsG0,chi0) 110 111 !Arguments ------------------------------------ 112 !scalars 113 logical,intent(in) :: is_metallic 114 integer,intent(in) :: ik_bz,npwepG0,nspinor 115 type(gsphere_t),intent(in) :: Gsph_epsG0 116 type(littlegroup_t),intent(in) :: Ltg_q 117 type(em1params_t),intent(in) :: Ep 118 !arrays 119 complex(gwpc),intent(in) :: rhotwg(npwepG0) 120 complex(dpc),intent(in) :: green_w(Ep%nomega) 121 complex(gwpc),intent(inout) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega) 122 123 !Local variables------------------------------- 124 !scalars 125 integer :: itim,io,isym,nthreads 126 integer :: isymop,nsymop 127 real(gwp) :: dr 128 complex(gwpc) :: dd 129 !character(len=500) :: msg 130 !arrays 131 integer :: Sm1_gmG0(Ep%npwe) 132 complex(gwpc),allocatable :: rhotwg_sym(:,:) 133 134 ! ************************************************************************* 135 136 ABI_UNUSED(nspinor) 137 138 nthreads = xomp_get_max_threads() 139 140 SELECT CASE (Ep%symchi) 141 CASE (0) 142 ! Do not use symmetries 143 144 ! note that single precision is faster (sometimes factor ~2). 145 ! Rely on MKL threads for OPENMP parallelization 146 147 do io=1,Ep%nomega 148 ! Check if green_w(io) is real (=> pure imaginary omega) 149 ! and that it is not a metal 150 ! then the corresponding chi0(io) is hermitian 151 if( ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic ) then 152 dr=green_w(io) 153 call xher('U',Ep%npwe,dr,rhotwg,1,chi0(:,:,io),Ep%npwe) 154 else 155 dd=green_w(io) 156 call xgerc(Ep%npwe,Ep%npwe,dd,rhotwg,1,rhotwg,1,chi0(:,:,io),Ep%npwe) 157 endif 158 end do 159 160 CASE (1) 161 ! Use symmetries to reconstruct the integrand in the BZ. 162 ! 163 ! Notes on the symmetrization of the oscillator matrix elements 164 ! If Sq = q then M_G^( Sk,q)= e^{-i(q+G).t} M_{ S^-1G} (k,q) 165 ! If -Sq = q then M_G^(-Sk,q)= e^{-i(q+G).t} M_{-S^-1G}^*(k,q) 166 ! 167 ! In case of an umklapp process 168 ! If Sq = q+G0 then M_G( Sk,q)= e^{-i(q+G).t} M_{ S^-1(G-G0} (k,q) 169 ! If -Sq = q+G0 then M_G(-Sk,q)= e^{-i(q+G).t} M_{-S^-1(G-G0)}^*(k,q) 170 ! 171 ! Ltg_q%igmG0(ig,itim,isym) contains the index of G-G0 where ISq=q+G0 172 ! Note that there is no need to take into account the phases due to q, 173 ! They cancel in the scalar product ==> phmGt(G,isym)=e^{-iG\cdot t} 174 ! 175 ! Mind the slicing of %rottbm1(npwepG0,timrev,nsym) and %phmGt(npwepG0,nsym) as 176 ! these arrays, usually, do not conform to rho_twg_sym(npw) ! 177 ! 178 ! Loop over symmetries of the space group and time-reversal. 179 nsymop = count(Ltg_q%wtksym(:,:,ik_bz)==1) 180 ABI_MALLOC(rhotwg_sym,(Ep%npwe,nsymop)) 181 isymop = 0 182 183 ! Prepare all the rhotwg at once to use BLAS level 3 routines 184 do isym=1,Ltg_q%nsym_sg 185 do itim=1,Ltg_q%timrev 186 if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then 187 ! This operation belongs to the little group and has to be used to reconstruct the BZ === 188 ! * In the following 3 lines mind the slicing (1:npwe) 189 ! TODO this is a hot-spot, should add a test on the umklapp 190 ! 191 !gmG0 => Ltg_q%igmG0(1:Ep%npwe,itim,isym) 192 Sm1_gmG0(1:Ep%npwe) = Gsph_epsG0%rottbm1( Ltg_q%igmG0(1:Ep%npwe,itim,isym), itim,isym) 193 194 isymop = isymop + 1 195 SELECT CASE (itim) 196 CASE (1) 197 rhotwg_sym(1:Ep%npwe,isymop) = rhotwg(Sm1_gmG0) * Gsph_epsG0%phmGt(1:Ep%npwe,isym) 198 CASE (2) 199 rhotwg_sym(1:Ep%npwe,isymop) = GWPC_CONJG(rhotwg(Sm1_gmG0))*Gsph_epsG0%phmGt(1:Ep%npwe,isym) 200 CASE DEFAULT 201 ABI_BUG(sjoin('Wrong itim:', itoa(itim))) 202 END SELECT 203 end if 204 end do 205 end do 206 207 ! Multiply rhotwg_sym by green_w(io) and accumulate in chi0(G,Gp,io) 208 ! note that single precision is faster (sometimes factor ~2). 209 ! Rely on MKL threads for OPENMP parallelization 210 do io=1,Ep%nomega 211 ! Check if green_w(io) is real (=> pure imaginary omega) 212 ! and that it is not a metal 213 ! then the corresponding chi0(io) is hermitian 214 if( ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic ) then 215 dr=green_w(io) 216 call xherk('U','N',Ep%npwe,nsymop,dr,rhotwg_sym,Ep%npwe,one_gw,chi0(:,:,io),Ep%npwe) 217 else 218 dd=green_w(io) 219 call xgemm('N','C',Ep%npwe,Ep%npwe,nsymop,dd,rhotwg_sym,Ep%npwe,rhotwg_sym,Ep%npwe,cone_gw,chi0(:,:,io),Ep%npwe) 220 endif 221 end do 222 223 ABI_FREE(rhotwg_sym) 224 225 CASE DEFAULT 226 ABI_BUG(sjoin('Wrong symchi:', itoa(Ep%symchi))) 227 END SELECT 228 229 end subroutine assemblychi0_sym
m_chi0tk/assemblychi0sf [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
assemblychi0sf
FUNCTION
Update the spectral function of the irreducible polarizability for the contribution of one pair of occupied-unoccupied states, for each frequency. If symchi==1, the symmetries of the little group of the external q-point are used to symmetrize the contribution in the full Brillouin zone. In this case, the routine computes: $ chi0(G1,G2,io)=chi0(G1,G2,io)+\sum_S (rhotwg(G1)*rhotwg^\dagger(G2))*\delta(w - trans) $ where S are the symmetries of the little group of the external q-point.
INPUTS
ik_bz=Index of the k-point in the BZ whose contribution has to be added to the spectral function of chi0 If symchi=1, the contribution is symmetrized. my_wl,my_wr=min and Max frequency index treated by this processor. npwe=Number of plane waves used to describe chi0. npwepG0=Maximum number of G vectors taking into account umklapp vectors. nomegasf=Number of frequencies for the spectral function. nspinor=Number of spinorial components. symchi=1 if symmetries are used, 0 otherwise rhotwg(npwepG0)=Oscillator matrix elements corresponding to an occupied-unoccupied pair of states. timrev=if 2, time reversal has to be used to obtain k_bz; 1 otherwise. Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors %ng=number of G vectors in the enlarged sphere, actually MUST be equal to the size of rhotwg %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion %phmGt(ng,nsym)=phase factors associated to non-simmorphic operations Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point. %timrev=2 it time-reversal is used, 1 otherwise %nsym_sg=Number of space group symmetries %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point %flag_umklp(timrev,nsym)= flag for umklapp processes if 1 that the particular operation (IS) requires a G_o to preserve Q, 0 otherwise %igmG0(npwepG0,timrev,nsym) index of G-G0 in the array gvec factocc=occupation factor=f_occ*(ockp-occk) (see cchi0.F90) wl,wr=Weights used to approximate the delta function.
OUTPUT
(see side effects)
SIDE EFFECTS
chi0sf(npwe,npwe,my_wl:my_wr)= updated spectral function.
NOTES
Umklapp processes are not yet implemented
SOURCE
1246 subroutine assemblychi0sf(ik_bz,symchi,Ltg_q,npwepG0,npwe,rhotwg,Gsph_epsG0,& 1247 & factocc,my_wl,iomegal,wl,my_wr,iomegar,wr,nomegasf,chi0sf) 1248 1249 !Arguments ------------------------------------ 1250 !scalars 1251 integer,intent(in) :: ik_bz,iomegal,iomegar,my_wl,my_wr,nomegasf,npwe,npwepG0 1252 integer,intent(in) :: symchi 1253 real(dp),intent(in) :: factocc,wl,wr 1254 type(gsphere_t),intent(in) :: Gsph_epsG0 1255 type(littlegroup_t),intent(in) :: Ltg_q 1256 !arrays 1257 complex(gwpc),intent(in) :: rhotwg(npwepG0) 1258 complex(gwpc),intent(inout) :: chi0sf(npwe,npwe,my_wl:my_wr) 1259 1260 !Local variables------------------------------- 1261 !scalars 1262 integer :: isym,itim,ig1,ig2 1263 complex(gwpc) :: num 1264 character(len=500) :: msg 1265 !arrays 1266 integer :: Sm1_gmG0(npwe) 1267 complex(gwpc) :: rhotwg_sym(npwe) 1268 1269 ! ************************************************************************* 1270 1271 if (iomegal < my_wl .or. iomegar > my_wr) then 1272 write(msg,'(3a,2(a,i0,a,i0,a))')ch10,& 1273 & ' Indices out of boundary ',ch10,& 1274 & ' my_wl = ',my_wl,' iomegal = ',iomegal,ch10,& 1275 & ' my_wr = ',my_wr,' iomegar = ',iomegar,ch10 1276 ABI_BUG(msg) 1277 end if 1278 1279 SELECT CASE (symchi) 1280 CASE (0) 1281 ! Do not use symmetries. 1282 1283 ! MG: This is the best I can do for this part. 1284 !$omp PARALLEL private(num) 1285 !$omp SECTIONS 1286 !$omp SECTION 1287 if (wl<huge(zero)*1.d-11) then !FIXME this is awful 1288 num=-wl*factocc 1289 call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,chi0sf(:,:,iomegal),npwe) 1290 end if 1291 1292 ! Last point, must accumulate left point but not the right one 1293 !$omp SECTION 1294 if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then 1295 num=-wr*factocc 1296 call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,chi0sf(:,:,iomegar),npwe) 1297 end if 1298 !$omp end SECTIONS 1299 !$omp end PARALLEL 1300 1301 CASE (1) 1302 ! Use symmetries to reconstruct oscillator matrix elements 1303 ! Notes on the symmetrization of the oscillator maxtri elements: 1304 ! 1305 ! If Sq=q then M_G^( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1G} (k,q) 1306 ! If -Sq=q then M_G^(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1G}^*(k,q) 1307 ! 1308 ! In case of an umklapp process 1309 ! If Sq=q+G_o then M_G( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1(G-G_o} (k,q) 1310 ! If -Sq=q+G_o then M_G(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1(G-G_o)}^*(k,q) 1311 ! 1312 ! Ltg_q%igmG0(ig,itim,isym) contains the index of G-G0 where ISq=q+G0 1313 ! Note that there is no need to take into account the phases due to q, 1314 ! They cancel in the scalar product ==> phmGt(G,isym)=e^{-iG\cdot t} 1315 ! 1316 ! Mind the slicing of %rottbm1(npwepG0,timrev,nsym) and %phgt(npwepG0,nsym) as 1317 ! these arrays, usually, do not conform to rho_twg_sym(npw) ! 1318 ! 1319 !ABI_MALLOC(rhotwg_sym,(npwe)) 1320 ! 1321 ! Loop over symmetries of the space group and time-reversal 1322 do isym=1,Ltg_q%nsym_sg 1323 do itim=1,Ltg_q%timrev 1324 1325 if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then 1326 ! This operation belongs to the little group and has to be used to reconstruct BZ. 1327 ! TODO this is a hot-spot, should add a test on the umklapp 1328 ! 1329 ! In these 3 lines mind the slicing (1:npwe) 1330 Sm1_gmG0(1:npwe)=Gsph_epsG0%rottbm1( Ltg_q%igmG0(1:npwe,itim,isym), itim,isym) 1331 1332 SELECT CASE (itim) 1333 CASE (1) 1334 rhotwg_sym(1:npwe)=rhotwg(Sm1_gmG0(1:npwe)) * Gsph_epsG0%phmGt(1:npwe,isym) 1335 CASE (2) 1336 rhotwg_sym(1:npwe)=conjg(rhotwg(Sm1_gmG0(1:npwe))) * Gsph_epsG0%phmGt(1:npwe,isym) 1337 CASE DEFAULT 1338 ABI_BUG(sjoin('Wrong value for itim:', itoa(itim))) 1339 END SELECT 1340 1341 #if 0 1342 !! MG: This is the best I can do, at present. 1343 !$omp PARALLEL private(num) 1344 !$omp SECTIONS 1345 1346 !$omp SECTION 1347 if (wl<huge(zero)*1.d-11) then 1348 num=-wl*factocc 1349 call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegal),npwe) 1350 end if 1351 !$omp SECTION 1352 ! 1353 ! Last point, must accumulate left point but not the right one 1354 if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then 1355 num=-wr*factocc 1356 call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegar),npwe) 1357 end if 1358 !$omp end SECTIONS 1359 !$omp end PARALLEL 1360 #else 1361 1362 if (wl<huge(zero)*1.d-11) then 1363 !call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegal),npwe) 1364 num=-wl*factocc 1365 !$omp parallel do 1366 do ig2=1,npwe 1367 do ig1=1,npwe 1368 chi0sf(ig1,ig2,iomegal) = chi0sf(ig1,ig2,iomegal) + num * rhotwg_sym(ig1) * conjg(rhotwg_sym(ig2)) 1369 end do 1370 end do 1371 end if 1372 1373 ! Last point, must accumulate left point but not the right one 1374 if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then 1375 !call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegar),npwe) 1376 num=-wr*factocc 1377 !$omp parallel do 1378 do ig2=1,npwe 1379 do ig1=1,npwe 1380 !call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegal),npwe) 1381 chi0sf(ig1,ig2,iomegar) = chi0sf(ig1,ig2,iomegar) + num * rhotwg_sym(ig1) * conjg(rhotwg_sym(ig2)) 1382 end do 1383 end do 1384 end if 1385 #endif 1386 end if !wtksym 1387 1388 end do !inv 1389 end do !isym 1390 !ABI_FREE(rhotwg_sym) 1391 1392 CASE DEFAULT 1393 ABI_BUG(sjoin('Wrong value for symchi:', itoa(symchi))) 1394 END SELECT 1395 1396 end subroutine assemblychi0sf
m_chi0tk/calc_kkweight [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
calc_kkweight
FUNCTION
Calculate frequency dependent weights needed to perform the Hilbert transform Subroutine needed to implement the calculation of the polarizability using the spectral representation as proposed in: PRB 74, 035101 (2006) [[cite:Shishkin2006]] and PRB 61, 7172 (2000) [[cite:Miyake2000]]
INPUTS
nsp=number of frequencies where the imaginary part of the polarizability is evaluated ne=number of frequencies for the polarizability (same as in epsilon^-1) omegasp(nsp)=real frequencies for the imaginary part of the polarizability omegae(ne)= imaginary frequencies for the polarizability delta=small imaginary part used to avoid poles, input variables
OUTPUT
kkweight(nsp,ne)=frequency dependent weights Eq A1 PRB 74, 035101 (2006) [[cite:Shishkin2006]]
SOURCE
1513 subroutine calc_kkweight(ne,omegae,nsp,omegasp,delta,omegamax,kkw) 1514 1515 !Arguments ------------------------------------ 1516 !scalars 1517 integer,intent(in) :: ne,nsp 1518 real(dp),intent(in) :: delta,omegamax 1519 !arrays 1520 real(dp),intent(in) :: omegasp(nsp) 1521 complex(dpc),intent(in) :: omegae(ne) 1522 complex(dpc),intent(out) :: kkw(nsp,ne) 1523 1524 !Local variables------------------------------- 1525 !scalars 1526 integer :: isp,je 1527 real(dp) :: eta,xx1,xx2,den1,den2 1528 complex(dpc) :: c1,c2,wt 1529 !************************************************************************ 1530 1531 DBG_ENTER("COLL") 1532 1533 kkw(:,:)=czero 1534 1535 do je=1,ne 1536 eta=delta 1537 wt=omegae(je) 1538 ! Not include shift at omega==0, what about metallic systems? 1539 if (abs(real(omegae(je)))<tol6 .and. abs(aimag(wt))<tol6) eta=tol12 1540 ! Not include shift along the imaginary axis 1541 if (abs(aimag(wt))>tol6) eta=zero 1542 do isp=1,nsp 1543 if (isp==1) then 1544 ! Skip negative point, should check that this would not lead to spurious effects 1545 c1=czero 1546 den1=one 1547 else 1548 xx1=omegasp(isp-1) 1549 xx2=omegasp(isp) 1550 den1= xx2-xx1 1551 c1= -(wt-xx1+j_dpc*eta)*log( (wt-xx2+j_dpc*eta)/(wt-xx1+j_dpc*eta) )& 1552 & +(wt+xx1-j_dpc*eta)*log( (wt+xx2-j_dpc*eta)/(wt+xx1-j_dpc*eta) ) 1553 c1= c1/den1 1554 end if 1555 xx1=omegasp(isp) 1556 if (isp==nsp) then 1557 ! Skip last point should check that this would not lead to spurious effects 1558 xx2=omegamax 1559 else 1560 xx2=omegasp(isp+1) 1561 end if 1562 den2=xx2-xx1 1563 c2= (wt-xx2+j_dpc*eta)*log( (wt-xx2+j_dpc*eta)/(wt-xx1+j_dpc*eta) )& 1564 & -(wt+xx2-j_dpc*eta)*log( (wt+xx2-j_dpc*eta)/(wt+xx1-j_dpc*eta) ) 1565 c2= c2/den2 1566 kkw(isp,je)= c1/den1 + c2/den2 1567 end do 1568 end do 1569 1570 DBG_EXIT("COLL") 1571 1572 end subroutine calc_kkweight
m_chi0tk/chi0_bbp_mask [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
chi0_bbp_mask
FUNCTION
INPUTS
OUTPUT
SOURCE
2338 subroutine chi0_bbp_mask(ikmq_ibz, ik_ibz, spin, spin_fact, use_tr, & 2339 gwcomp, spmeth, chi_nband, mband, ebands, bbp_mask) 2340 2341 !Arguments ------------------------------------ 2342 !scalars 2343 integer,intent(in) :: spin,ik_ibz,ikmq_ibz,mband,gwcomp, spmeth, chi_nband 2344 real(dp),intent(in) :: spin_fact 2345 logical,intent(in) :: use_tr 2346 type(ebands_t),target,intent(in) :: ebands 2347 !arrays 2348 logical,intent(out) :: bbp_mask(mband, mband) 2349 2350 !Local variables------------------------------- 2351 !scalars 2352 integer :: ib1, ib2 2353 real(dp) :: deltaeGW_b1kmq_b2k,deltaf_b1kmq_b2k,e_b1_kmq,f_b1_kmq 2354 !arrays 2355 real(dp), ABI_CONTIGUOUS pointer :: qp_eig(:,:,:),qp_occ(:,:,:) 2356 2357 !************************************************************************ 2358 2359 qp_eig => ebands%eig; qp_occ => ebands%occ 2360 bbp_mask = .FALSE. 2361 2362 select case (gwcomp) 2363 case (0) 2364 ! Loop over "conduction" states. 2365 do ib1=1,chi_nband 2366 e_b1_kmq = qp_eig(ib1, ikmq_ibz, spin) 2367 f_b1_kmq = qp_occ(ib1, ikmq_ibz, spin) 2368 2369 ! Loop over "valence" states. 2370 do ib2=1,chi_nband 2371 deltaf_b1kmq_b2k = spin_fact * (f_b1_kmq - qp_occ(ib2,ik_ibz,spin)) 2372 deltaeGW_b1kmq_b2k = e_b1_kmq - qp_eig(ib2,ik_ibz,spin) 2373 2374 select case (spmeth) 2375 case (0) 2376 ! Standard Adler-Wiser expression. 2377 if (ABS(deltaf_b1kmq_b2k) >= GW_TOL_DOCC) then 2378 bbp_mask(ib1, ib2) = .TRUE. 2379 if (use_tr .and. ib1 < ib2) bbp_mask(ib1,ib2) = .FALSE. ! GAIN a factor ~2 thanks to time-reversal. 2380 end if 2381 2382 case (1,2) 2383 ! Spectral method, WARNING time-reversal here is always assumed! 2384 if (ABS(deltaf_b1kmq_b2k) >= GW_TOL_DOCC) then 2385 bbp_mask(ib1,ib2)=.TRUE. 2386 if (deltaeGW_b1kmq_b2k<zero) bbp_mask(ib1,ib2)=.FALSE. ! Only positive frequencies are needed for the Hilbert transform. 2387 !$if (use_tr .and. ib1<ib2) bbp_mask(ib1,ib2)=.FALSE. ! GAIN a factor ~2 thanks to time-reversal. 2388 end if 2389 2390 case default 2391 ABI_ERROR(sjoin("Wrong value for spmeth:", itoa(spmeth))) 2392 end select 2393 !write(std_out,*) "bbp_mask(ib1,ib2)",bbp_mask(ib1,ib2) 2394 end do !ib2 2395 end do !ib1 2396 2397 case (1) 2398 ! Extrapolar technique 2399 ABI_CHECK(spmeth == 0, "Hilbert transform and extrapolar method are not compatible") 2400 2401 ! Loop over "conduction" states. 2402 do ib1=1,chi_nband 2403 e_b1_kmq=qp_eig(ib1,ikmq_ibz,spin) 2404 f_b1_kmq= qp_occ(ib1,ikmq_ibz,spin) 2405 2406 ! Loop over "valence" states. 2407 do ib2=1,chi_nband 2408 deltaf_b1kmq_b2k = spin_fact*(f_b1_kmq-qp_occ(ib2,ik_ibz,spin)) 2409 deltaeGW_b1kmq_b2k= e_b1_kmq-qp_eig(ib2,ik_ibz,spin) 2410 2411 ! When the completeness correction is used, 2412 ! we need to also consider transitions with vanishing deltaf 2413 ! Rangel: This is to compute chi in metals correctly with the extrapolar method. 2414 bbp_mask(ib1,ib2)=.TRUE. 2415 !if (qp_occ(ib2,ik_ibz,is) < GW_TOL_DOCC) CYCLE 2416 if (qp_occ(ib2,ik_ibz,spin) < GW_TOL_DOCC .and. (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC .or. ib1<ib2)) then 2417 bbp_mask(ib1,ib2)=.FALSE. 2418 end if 2419 end do 2420 end do 2421 2422 case default 2423 ABI_ERROR(sjoin("Wrong value of gwcomp:", itoa(gwcomp))) 2424 end select 2425 2426 end subroutine chi0_bbp_mask
m_chi0tk/completechi0_deltapart [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
completechi0_deltapart
FUNCTION
Apply the delta part of the completeness correction to chi0
INPUTS
ik_bz=Index of the k-point in the full BZ whose contribution has to be added and symmetrized. qzero=.TRUE. is long wave-length limit. symchi=1 if we are summing over IBZ_q and symmetrization has to be performed. npwe=Number of G vectors in chi0. npwvec=MAX number of G. nomega=Number of frequencies. nspinor=Number of spinorial components. nfftot=Total Number of points in the FFT ngfft(18)=Info on the FFT. igfft0(npwvec)=Index of each G in the FFT array. Gsph_FFT=<gsphere_t>=Info on the largest G-sphere contained in the FFT box used for wavefunctions. Ltg_q=<littlegroup_t>= Structure gathering information on the little group of the external q. green_enhigh_w=Approximated frequency dependent part of the Green function entering equation (TODO put reference) wfwfg=Fourier components of u_{kb1}.u_{kb2}
OUTPUT
See SIDES EFFECTS SIDES EFFECTS chi0(npwe,npwe,nomega)= In input chi0 calculated so far, In output the "delta part" of the completeness correction is added.
SOURCE
1911 subroutine completechi0_deltapart(ik_bz,qzero,symchi,npwe,npwvec,nomega,nspinor,& 1912 & nfftot,ngfft,igfft0,Gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0) 1913 1914 !Arguments ------------------------------------ 1915 !scalars 1916 integer,intent(in) :: ik_bz,nfftot,nomega,npwe,npwvec,nspinor,symchi 1917 logical,intent(in) :: qzero 1918 type(gsphere_t),intent(in) :: Gsph_FFT 1919 type(littlegroup_t),intent(in) :: Ltg_q 1920 !arrays 1921 integer,intent(in) :: igfft0(npwvec),ngfft(18) 1922 complex(dpc),intent(in) :: green_enhigh_w(nomega) 1923 complex(gwpc),intent(in) :: wfwfg(nfftot*nspinor**2) 1924 complex(gwpc),intent(inout) :: chi0(npwe,npwe,nomega) 1925 1926 !Local variables ------------------------------ 1927 !scalars 1928 integer,save :: enough=0 1929 integer :: iSm1_g1mg2,iSm1_g1mg2_fft,ig,gmg_sph,gmg_fft 1930 integer :: igp,igstart,isym,itim,outofbox_wfn 1931 complex(gwpc) :: phmGt 1932 !character(len=500) :: msg 1933 1934 !************************************************************************ 1935 1936 igstart=1; if (qzero) igstart=2 1937 outofbox_wfn=0 1938 1939 SELECT CASE (symchi) 1940 1941 CASE (0) ! Do not use symmetries. 1942 ! MG: One has to make sure G1-G2 is still in the FFT mesh for each G1 and G2 in chi0 (not always true) 1943 ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic! 1944 do igp=igstart,npwe 1945 do ig=igstart,npwe 1946 gmg_fft = Gsph_FFT%gmg_fftidx(ig,igp,ngfft) 1947 if (gmg_fft==0) then 1948 outofbox_wfn=outofbox_wfn+1; CYCLE 1949 end if 1950 chi0(ig,igp,:) = chi0(ig,igp,:) + wfwfg(gmg_fft)*green_enhigh_w(:) 1951 end do 1952 end do 1953 1954 CASE (1) 1955 ! Symmetrize the integrand in the full BZ. 1956 ! * <Sk b|e^{-i(G1-G2}.r}|b Sk> = e^{-i(G1-G2).\tau} <k b|e^{-i(S^{-1}(G1-G2).r)|b k> 1957 ! * green_enhigh_w in invariant under symmetry 1958 ! * We symmetrize using the operations of the little group of q since this routine 1959 ! is called inside a sum over IBZ_q, it would be possible to symmetrize 1960 ! this term by just summing over the IBZ and rotating the matrix elements. 1961 ! * Time-reversal does not lead to a complex conjugated since bra and ket are the same. 1962 ! 1963 do igp=igstart,npwe 1964 do ig=igstart,npwe 1965 1966 ! Get the index of G1-G2. 1967 gmg_sph = Gsph_FFT%gmg_idx(ig,igp) 1968 if (gmg_sph==0) then 1969 outofbox_wfn=outofbox_wfn+1; CYCLE 1970 end if 1971 1972 do itim=1,Ltg_q%timrev 1973 do isym=1,Ltg_q%nsym_sg 1974 if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then 1975 ! * This operation belongs to the little group and has to be used to reconstruct the BZ. 1976 ! * Time-reversal in not used to rotate (G1-G2) see comment above. 1977 phmGt = Gsph_FFT%phmGt (gmg_sph,isym) 1978 iSm1_g1mg2 = Gsph_FFT%rottbm1(gmg_sph,1,isym) 1979 iSm1_g1mg2_fft = igfft0(iSm1_g1mg2) 1980 1981 chi0(ig,igp,:) = chi0(ig,igp,:) + phmGt*wfwfg(iSm1_g1mg2_fft)*green_enhigh_w(:) 1982 end if 1983 end do !isym 1984 end do !itim 1985 1986 end do !igp 1987 end do !ig 1988 1989 CASE DEFAULT 1990 ABI_BUG("Wrong value of symchi") 1991 END SELECT 1992 1993 if (outofbox_wfn/=0) then 1994 enough=enough+1 1995 if (enough<=50) then 1996 ABI_WARNING(sjoin(' Number of G1-G2 pairs outside the G-sphere for Wfns: ', itoa(outofbox_wfn))) 1997 if (enough==50) then 1998 call wrtout(std_out,' ========== Stop writing Warnings ==========') 1999 end if 2000 end if 2001 end if 2002 2003 end subroutine completechi0_deltapart
m_chi0tk/hilbert_transform [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
hilbert_transform
FUNCTION
Compute the hilbert transform.
INPUTS
nomegasf=number of points for the imaginary part of $\chi0(q,\omega)$ nomega=number of frequencies in $\chi0(q,\omega)$. max_rest,min_res=max and min resonant transition energy (for this q-point) my_max_rest,my_min_rest=max and min resonant transition energy treated by this processor
OUTPUT
SOURCE
1746 subroutine hilbert_transform(npwe,nomega,nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,spmeth) 1747 1748 !Arguments ------------------------------------ 1749 !scalars 1750 integer,intent(in) :: spmeth,nomega,nomegasf,my_wl,my_wr,npwe 1751 !arrays 1752 complex(dpc),intent(in) :: kkweight(nomegasf,nomega) 1753 complex(gwpc),intent(inout) :: sf_chi0(npwe,npwe,my_wl:my_wr) 1754 complex(gwpc),intent(inout) :: chi0(npwe,npwe,nomega) 1755 1756 !Local variables------------------------------- 1757 !scalars 1758 integer :: ig2,my_nwp 1759 character(len=500) :: msg 1760 !arrays 1761 complex(gwpc),allocatable :: A_g1wp(:,:),H_int(:,:),my_kkweight(:,:) 1762 1763 !************************************************************************ 1764 1765 #ifdef HAVE_OPENMP 1766 write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform (with OpenMP) using method ',spmeth,' It might take some time...' 1767 #else 1768 write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform using method ',spmeth,' It might take some time...' 1769 #endif 1770 call wrtout(std_out, msg, do_flush=.True.) 1771 1772 my_nwp = my_wr - my_wl +1 1773 1774 !$omp parallel private(my_kkweight, A_g1wp, H_int, ig2) 1775 ABI_MALLOC(my_kkweight, (my_wl:my_wr,nomega)) 1776 my_kkweight = kkweight(my_wl:my_wr,:) 1777 1778 ABI_MALLOC(A_g1wp, (npwe, my_nwp)) 1779 ABI_MALLOC(H_int, (npwe, nomega)) 1780 1781 !$omp do 1782 do ig2=1,npwe 1783 A_g1wp = sf_chi0(:,ig2,:) 1784 1785 ! Compute H_int = MATMUL(A_g1wp,my_kkweight) 1786 call XGEMM('N','N',npwe,nomega,my_nwp,cone_gw,A_g1wp,npwe,my_kkweight,my_nwp,czero_gw,H_int,npwe) 1787 chi0(:,ig2,:) = H_int 1788 end do 1789 1790 ABI_FREE(my_kkweight) 1791 ABI_FREE(A_g1wp) 1792 ABI_FREE(H_int) 1793 !$omp end parallel 1794 1795 end subroutine hilbert_transform
m_chi0tk/hilbert_transform_headwings [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
hilbert_transform_headwings
FUNCTION
Compute the hilbert transform the heads and wings of the polarizability.
INPUTS
nomegasf=number of points for the imaginary part of $\chi0(q,\omega)$ nomega=number of frequencies in $\chi0(q,\omega)$. max_rest,min_res=max and min resonant transition energy (for this q-point) my_max_rest,my_min_rest=max and min resonant transition energy treated by this processor
OUTPUT
SOURCE
1817 subroutine hilbert_transform_headwings(npwe,nomega,nomegasf,my_wl,my_wr,kkweight, & 1818 & sf_lwing,sf_uwing,sf_head,chi0_lwing,chi0_uwing,chi0_head,spmeth) 1819 1820 !Arguments ------------------------------------ 1821 !scalars 1822 integer,intent(in) :: spmeth,nomega,nomegasf,my_wl,my_wr,npwe 1823 !arrays 1824 complex(dpc),intent(in) :: kkweight(nomegasf,nomega) 1825 complex(dpc),intent(inout) :: sf_lwing(npwe,my_wl:my_wr,3) 1826 complex(dpc),intent(inout) :: sf_uwing(npwe,my_wl:my_wr,3) 1827 complex(dpc),intent(inout) :: sf_head(3,3,my_wl:my_wr) 1828 complex(dpc),intent(inout) :: chi0_lwing(npwe,nomega,3) 1829 complex(dpc),intent(inout) :: chi0_uwing(npwe,nomega,3) 1830 complex(dpc),intent(inout) :: chi0_head(3,3,nomega) 1831 1832 !Local variables------------------------------- 1833 !scalars 1834 integer :: ig1,idir,io,iw 1835 complex(dpc) :: kkw 1836 character(len=500) :: msg 1837 !************************************************************************ 1838 1839 #ifdef HAVE_OPENMP 1840 write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform (with OpenMP) using method ',spmeth,' It might take some time...' 1841 #else 1842 write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform using method ',spmeth,' It might take some time...' 1843 #endif 1844 call wrtout(std_out,msg, do_flush=.True.) 1845 1846 ! Hilbert transform of the head. 1847 do io=1,nomega 1848 chi0_head(1,1,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(1,1,my_wl:my_wr)) 1849 chi0_head(2,1,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(2,1,my_wl:my_wr)) 1850 chi0_head(3,1,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(3,1,my_wl:my_wr)) 1851 chi0_head(1,2,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(1,2,my_wl:my_wr)) 1852 chi0_head(2,2,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(2,2,my_wl:my_wr)) 1853 chi0_head(3,2,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(3,2,my_wl:my_wr)) 1854 chi0_head(1,3,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(1,3,my_wl:my_wr)) 1855 chi0_head(2,3,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(2,3,my_wl:my_wr)) 1856 chi0_head(3,3,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(3,3,my_wl:my_wr)) 1857 end do 1858 1859 ! Hilbert transform for wings. 1860 ! Partial contributions to chi0 will be summed afterwards. 1861 !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(kkw) 1862 do idir=1,3 1863 do io=1,nomega 1864 do iw=my_wl,my_wr 1865 kkw = kkweight(iw,io) 1866 do ig1=1,npwe 1867 chi0_lwing(ig1,io,idir) = chi0_lwing(ig1,io,idir) + kkw*sf_lwing(ig1,iw,idir) 1868 chi0_uwing(ig1,io,idir) = chi0_uwing(ig1,io,idir) + kkw*sf_uwing(ig1,iw,idir) 1869 end do 1870 end do 1871 end do 1872 end do !idir 1873 1874 end subroutine hilbert_transform_headwings
m_chi0tk/make_transitions [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
make_transitions
FUNCTION
Calculate transition energies entering the espression for the irreducible polarizability.
INPUTS
nsspol=1 for spin unpolarized, 2 for spin polarized calculations nbnds=total number of bands kmesh<kmesh_t>=datatype gathering info on the k-mesh: | %nbz=number of k-points in the full BZ | %nibz=number of k-points in the IBZ | %tab(nkbz)=table giving for each k-point in the BZ, the corresponding irreducible point in the IBZ array | %bz(3,nkbz)=reduced coordinated of k-points TOL_DELTA_OCC=tolerance on the difference of the occupation numbers gw_energy(nbnds,kmesh%nkibz,nsppol)=quasi-particle energies energies occ(nbnds,kmesh%nkibz,nsppol)=occupation numbers chi0alg=integer defining the method used to calculate chi0 0 ==> calculate chi0 using the Adler-Wiser expression 1 ==> use spectral method timrev=if 2, time-reversal symmetry is considered; 1 otherwise
OUTPUT
my_max_rest,my_min_rest=Maximum and minimum resonant (posite) transition energy. max_rest,min_rest=Maximun and minimum resonant (posite) transition energy treated by this node.
SOURCE
2193 subroutine make_transitions(Wfd,chi0alg,nbnds,nbvw,nsppol,symchi,timrev,TOL_DELTA_OCC,& 2194 max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,gw_energy,occ,qpoint,bbp_ks_distrb) 2195 2196 !Arguments ------------------------------------ 2197 !scalars 2198 integer,intent(in) :: chi0alg,nbnds,nbvw,nsppol,symchi,timrev 2199 real(dp),intent(in) :: TOL_DELTA_OCC 2200 real(dp),intent(out) :: max_rest,min_rest, my_max_rest,my_min_rest 2201 type(kmesh_t),intent(in) :: Kmesh 2202 type(littlegroup_t),intent(in) :: Ltg_q 2203 type(wfdgw_t),intent(in) :: Wfd 2204 !arrays 2205 real(dp),intent(in) :: gw_energy(nbnds,Kmesh%nibz,nsppol) 2206 real(dp),intent(in) :: occ(nbnds,Kmesh%nibz,nsppol),qpoint(3) 2207 integer,intent(in) :: bbp_ks_distrb(Wfd%mband,Wfd%mband,Kmesh%nbz,Wfd%nsppol) 2208 2209 !Local variables------------------------------- 2210 !scalars 2211 integer :: ib1,ib2,ii,ik_bz,ik_ibz,ikmq_bz,ikmq_ibz,is,nt,ntrans,my_ntrans,iloop 2212 real(dp) :: delta_ene,delta_occ,spin_fact 2213 character(len=500) :: msg 2214 !arrays 2215 integer :: G0(3) 2216 real(dp) :: kmq(3) 2217 2218 !************************************************************************ 2219 2220 DBG_ENTER("COLL") 2221 2222 if (chi0alg < 0 .or. chi0alg >= 2) then 2223 ABI_BUG(sjoin('chi0alg:', itoa(chi0alg),' not allowed')) 2224 end if 2225 if (timrev /= 1 .and. timrev /= 2) then 2226 ABI_BUG(sjoin('timrev:', itoa(timrev),' not allowed')) 2227 end if 2228 2229 ABI_UNUSED(nbvw) 2230 ! 2231 ! In the first loop, we calculate total number of transitions for this q-point 2232 ! as well the min and max transition without taking into account distribution of bands. 2233 ! In the second iteration, we calculate the min and Max transition treated by this processor. 2234 ! 2235 spin_fact =half; if (nsppol == 2) spin_fact = one 2236 2237 my_max_rest=smallest_real; my_min_rest=greatest_real 2238 max_rest=smallest_real; min_rest=greatest_real 2239 2240 do iloop=1,2 2241 nt=0 2242 do ik_bz=1,Kmesh%nbz 2243 ik_ibz=Kmesh%tab(ik_bz) 2244 kmq(:)=Kmesh%bz(:,ik_bz)-qpoint(:) 2245 2246 if (symchi == 1) then 2247 if (Ltg_q%ibzq(ik_bz) /= 1) cycle ! This point does not belong to the IBZ defined by the little group 2248 end if 2249 2250 ! Find kp=k-q-G0 and also G0 where kp is in the first BZ 2251 if (.not. kmesh%has_BZ_item(kmq,ikmq_bz,g0)) then 2252 ! Stop as the weight 1.0/nkbz is wrong. 2253 write(msg,'(4a,2(2a,3f12.6),2a)')ch10,& 2254 ' make_transitions : ERROR - ',ch10,& 2255 ' kp = k-q-G0 not found in the BZ mesh',ch10,& 2256 ' k = ',(Kmesh%bz(ii,ik_bz),ii=1,3),ch10,& 2257 ' k-q = ',(kmq(ii),ii=1,3),ch10,& 2258 ' weight in cchi0/cchi0q is wrong ' 2259 ABI_ERROR(msg) 2260 end if 2261 2262 ikmq_ibz = Kmesh%tab(ikmq_bz) 2263 do is=1,nsppol 2264 do ib1=1,nbnds 2265 do ib2=1,nbnds 2266 2267 if (iloop == 2) then 2268 if (bbp_ks_distrb(ib1,ib2,ik_bz,is)/=Wfd%my_rank) cycle 2269 end if 2270 2271 if (timrev == 2 .and. ib1 < ib2) cycle ! Thanks to time-reversal we gain a factor ~2. 2272 2273 delta_occ = spin_fact * (occ(ib1,ikmq_ibz,is) - occ(ib2,ik_ibz,is)) 2274 delta_ene = gw_energy(ib1,ikmq_ibz,is) - gw_energy(ib2,ik_ibz,is) 2275 2276 if (chi0alg == 0) then 2277 ! Adler-Wiser expression. Skip only if factor due to occupation number is smaller than TOL_DELTA_OCC 2278 if (abs(delta_occ) < abs(TOL_DELTA_OCC)) cycle 2279 else if (chi0alg==1) then 2280 ! Spectral method with time-reversal, only resonant transitions 2281 ! This has to be changed to include spectral method without time-reversal 2282 if (delta_ene < -abs(TOL_DELTA_OCC) .or. abs(delta_occ) < abs(TOL_DELTA_OCC)) cycle 2283 end if 2284 2285 ! We have a new transition 2286 nt=nt+1 2287 2288 if (iloop==1) then 2289 max_rest=MAX(max_rest,zero,delta_ene) 2290 if (delta_ene>=-tol6) min_rest=MIN(min_rest,delta_ene) 2291 end if 2292 if (iloop==2) then 2293 my_max_rest=MAX(my_max_rest,zero,delta_ene) 2294 if (delta_ene>=-tol6) my_min_rest=MIN(my_min_rest,delta_ene) 2295 end if 2296 2297 end do 2298 end do 2299 end do 2300 end do 2301 if (iloop==1) ntrans=nt 2302 if (iloop==2) my_ntrans=nt 2303 end do !iloop 2304 2305 write(msg,'(2a,i9,2a,f8.3,3a,f8.3,a)')ch10,& 2306 ' Total number of transitions = ',ntrans,ch10,& 2307 ' min resonant = ',min_rest*Ha_eV,' [eV] ',ch10,& 2308 ' Max resonant = ',max_rest*Ha_eV,' [eV] ' 2309 call wrtout(std_out, msg) 2310 2311 if (Wfd%nproc/=1) then 2312 write(msg,'(2a,i9,2a,f8.3,3a,f8.3,a)')ch10,& 2313 ' Total number of transitions for this processor= ',my_ntrans,ch10,& 2314 ' min resonant = ',my_min_rest*Ha_eV,' [eV] ',ch10,& 2315 ' Max resonant = ',my_max_rest*Ha_eV,' [eV] ' 2316 call wrtout(std_out, msg) 2317 end if 2318 2319 DBG_EXIT("COLL") 2320 2321 end subroutine make_transitions
m_chi0tk/mkrhotwg_sigma [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
mkrhotwg_sigma
FUNCTION
Helper function used to calculate selected linear combination of the oscillator matrix elements in the case of non-collinear magnetism.
INPUTS
ii=Index selecting the particolar combination of spin components. npw=Number of plane-waves in the oscillators. nspinor=Number of spinorial components. rhotwg(npw*nspinor**2)=OScillator matrix elements.
OUTPUT
rhotwg_I(npw)=Required linear combination of the oscillator matrix elements.
SOURCE
253 subroutine mkrhotwg_sigma(ii,nspinor,npw,rhotwg,rhotwg_I) 254 255 !Arguments ------------------------------------ 256 !scalars 257 integer,intent(in) :: ii,npw,nspinor 258 !arrays 259 complex(gwpc),intent(in) :: rhotwg(npw*nspinor**2) 260 complex(gwpc),intent(out) :: rhotwg_I(npw) 261 262 ! ************************************************************************* 263 264 SELECT CASE (ii) 265 CASE (1) 266 ! $ M_0 = M_{\up,\up} + M_{\down,\down} $ 267 rhotwg_I(:) = rhotwg(1:npw) + rhotwg(npw+1:2*npw) 268 CASE (2) 269 ! $ M_z = M_{\up,\up} - M_{\down,\down} $ 270 rhotwg_I(:) = rhotwg(1:npw) - rhotwg(npw+1:2*npw) 271 CASE (3) 272 ! $ M_x = M_{\up,\down} + M_{\down,\up} $ 273 rhotwg_I(:) = ( rhotwg(2*npw+1:3*npw) + rhotwg(3*npw+1:4*npw) ) 274 CASE (4) 275 ! $ M_y = i * (M_{\up,\down} -M_{\down,\up}) $ 276 rhotwg_I(:) = (rhotwg(2*npw+1:3*npw) - rhotwg(3*npw+1:4*npw) )*j_gw 277 CASE DEFAULT 278 ABI_BUG(sjoin('Wrong ii value:', itoa(ii))) 279 END SELECT 280 281 end subroutine mkrhotwg_sigma
m_chi0tk/output_chi0sumrule [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
output_chi0sumrule
FUNCTION
Calculate and output the value of the sum rule for the non-interacting polarizability chi0
INPUTS
OUTPUT
(for writing routines, no output) otherwise, should be described
SOURCE
2024 subroutine output_chi0sumrule(qeq0,iq,npwe,omegaplasma,chi0sumrule,epsm1_w0,vc_sqrt) 2025 2026 !Arguments ------------------------------------ 2027 !scalars 2028 integer,intent(in) :: iq,npwe 2029 real(dp),intent(in) :: omegaplasma 2030 logical,intent(in) :: qeq0 2031 !arrays 2032 real(dp),intent(inout) :: chi0sumrule(npwe) 2033 complex(gwpc),intent(in) :: epsm1_w0(npwe,npwe),vc_sqrt(npwe) 2034 2035 !Local variables ------------------------------ 2036 !scalars 2037 integer :: ig,igstart 2038 real(dp) :: average,norm 2039 character(len=500) :: msg 2040 2041 !************************************************************************ 2042 2043 igstart=1; if (qeq0) igstart=2 2044 ! 2045 ! The sumrule reads: 2046 ! $ \int d\omega \omega v * Im[ \chi_0(\omega) ] = \pi/2 * w_p^2 $. 2047 chi0sumrule(igstart:npwe) = chi0sumrule(igstart:npwe) * vc_sqrt(igstart:npwe)**2 2048 ! 2049 ! Calculate a weighted average of the fulfilment of the sumrule on epsilon 2050 ! The weight is given according to the significance of each q+G in the 2051 ! subsequent GW calculation: It is proportional to v * (epsm1 -1 ) 2052 average = zero; norm = zero 2053 do ig=igstart,npwe 2054 average = average + chi0sumrule(ig) * real( vc_sqrt(ig)**2 * (epsm1_w0(ig,ig) - 1.0_dp ) ) 2055 norm = norm + real( vc_sqrt(ig)**2 * (epsm1_w0(ig,ig) - 1.0_dp ) ) 2056 !average = average + chi0sumrule(ig) * real( (epsm1_w0(ig,ig) - 1.0_dp ) ) 2057 !norm = norm + real( (epsm1_w0(ig,ig) - 1.0_dp ) ) 2058 !write(203,'(i4,8(2x,e12.6))') ig,1.0_dp/vc_sqrt(ig),chi0sumrule(ig)/ (0.5d0*omegaplasma**2*pi) 2059 end do 2060 2061 if (abs(norm)>tol8) then 2062 write(msg,'(1x,a,i4,a,f10.2,2x,a)')& 2063 ' Average fulfillment of the sum rule on Im[epsilon] for q-point ',& 2064 iq,' :',average/norm/(0.5_dp*omegaplasma**2*pi)*100.0_dp,'[%]' 2065 call wrtout([std_out, ab_out], msg) 2066 end if 2067 2068 end subroutine output_chi0sumrule
m_chi0tk/setup_spectral [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
setup_spectral
FUNCTION
Calculation of \chi_o based on the spectral method as proposed in PRB 74, 035101 (2006) [[cite:Shishkin2006]] and PRB 61, 7172 (2000) [[cite:Miyake2000]]. Setup of the real frequency mesh for $\Im\chi_o$ and of the frequency-dependent weights for Hilbert transform. Note that CPU time does not depend dramatically on nomegasf unlike memory. spmeth defines the approximant for the delta function: ==1 : use Triangular approximant (Kresse method) ==2 : use Gaussian method, requiring smearing (Miyake method)
INPUTS
nomegasf=number of points for the imaginary part of $\chi0(q,\omega)$ nomega=number of frequencies in $\chi0(q,\omega)$. max_rest,min_res=max and min resonant transition energy (for this q-point) my_max_rest,my_min_rest=max and min resonant transition energy treated by this processor method=integer flag defining the type of frequency mesh used for $\Im chi0$ | 0 for a linear mesh | 1 for a mesh densified around omegaplasma omegaplasma=frequency around which the mesh is densifies (usually Drude plasma frequency) used only in case of method==1 zcut=small imaginary shift to avoid pole in chi0
OUTPUT
kkweight(nomegasf,nomega)=Frequency dependent weight for Hilber transform. omegasf(nomegasf+1)=frequencies for imaginary part.
SOURCE
1608 subroutine setup_spectral(nomega,omega,nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,& 1609 & method,zcut,omegaplasma,my_wl,my_wr,kkweight) 1610 1611 !Arguments ------------------------------------ 1612 !scalars 1613 integer,intent(in) :: method,nomega,nomegasf 1614 integer,intent(out) :: my_wl,my_wr 1615 real(dp),intent(in) :: max_rest,min_rest,omegaplasma,zcut 1616 real(dp),intent(in) :: my_max_rest,my_min_rest 1617 !arrays 1618 real(dp),intent(out) :: omegasf(nomegasf) 1619 complex(dpc),intent(in) :: omega(nomega) 1620 complex(dpc),intent(out) :: kkweight(nomegasf,nomega) 1621 1622 !Local variables------------------------------- 1623 !scalars 1624 integer :: io,ii 1625 real(dp) :: nu_min,nu_max,nu1,nu2,dd,domegasf,wp,deltat 1626 character(len=500) :: msg 1627 !arrays 1628 integer,allocatable :: insort(:) 1629 !************************************************************************ 1630 1631 ! The mesh must enclose the entire range of transitions. 1632 dd=(max_rest-min_rest)/(nomegasf-1) 1633 domegasf=(max_rest-min_rest+2*dd)/(nomegasf-1) 1634 1635 write(msg,'(4a,f8.3,3a,i5,2a,f8.5,a)')ch10,& 1636 ' === Info on the real frequency mesh for spectral method === ',ch10,& 1637 ' maximum frequency = ',max_rest*Ha_eV,' [eV]',ch10,& 1638 ' nomegasf = ',nomegasf,ch10,& 1639 ' domegasf = ',domegasf*Ha_eV,' [eV]' 1640 call wrtout(std_out,msg) 1641 1642 if (min_rest<tol6) then 1643 ABI_WARNING("System seems to be metallic") 1644 end if 1645 1646 ! ====================================================== 1647 ! === Setup of the w-mesh for the spectral function ==== 1648 ! ====================================================== 1649 SELECT CASE (method) 1650 CASE (0) 1651 ! Linear mesh. 1652 call wrtout(std_out, ' Using linear mesh for Im chi0') 1653 do io=1,nomegasf 1654 omegasf(io)=(io-1)*domegasf+min_rest-dd 1655 end do 1656 1657 CASE (1) 1658 ! Non-homogeneous mesh densified around omega_plasma, do not improve results === 1659 ! WARNING_ this part has to be checked since I modified omegasf 1660 write(msg,'(a,f7.4,a)')' Using mesh densified around ',omegaplasma*Ha_eV,' [eV] ' 1661 call wrtout(std_out, msg) 1662 wp=omegaplasma ; deltat=max_rest-min_rest 1663 nu_min=zero 1664 if (deltat<wp ) then 1665 nu_max = wp/sqrt2 * ATAN(sqrt2*deltat*wp/(-deltat**2+wp**2)) 1666 else 1667 nu_max = wp/sqrt2 * ( ATAN(sqrt2*deltat*wp/(-deltat**2+wp**2)) + pi) 1668 end if 1669 domegasf=(nu_max-nu_min)/(nomegasf+1) 1670 !write(std_out,*) -(wp/sqrt2) * atan(sqrt2*deltat*wp/(deltat**2-wp**2)) 1671 omegasf(1)=zero ; omegasf(nomegasf+1)=deltat 1672 ii=0 1673 do io=2,nomegasf 1674 nu1=domegasf*(io-1) ; nu2=TAN(-sqrt2*nu1/wp) 1675 if (nu2<0) then 1676 omegasf(io) = wp * (one - SQRT(1+2*nu2**2))/(sqrt2*nu2) 1677 else 1678 omegasf(io) = wp * (one + SQRT(1+2*nu2**2))/(sqrt2*nu2) 1679 end if 1680 if (omegasf(io)> deltat ) then 1681 omegasf(io)= deltat-0.1*ii 1682 ii=ii+1 1683 end if 1684 ! write(102,'(i4,2x,3(f9.4,2x))')io,nu1,nu2,ep%omegasf(io)*Ha_eV 1685 end do 1686 1687 ! Reorder frequencies in ascending order 1688 ABI_MALLOC(insort,(nomegasf+1)) 1689 insort(:)=(/ (io,io=1,nomegasf+1) /) 1690 call sort_dp(nomegasf+1,omegasf,insort,tol14) 1691 ABI_FREE(insort) 1692 1693 CASE DEFAULT 1694 ABI_BUG(sjoin('Wrong value for method:', itoa(method))) 1695 END SELECT 1696 !write(std_out,*)omegasf(1)*Ha_eV,omegasf(nomegasf)*Ha_eV 1697 1698 ! Find min and max index in omegasf treated by this processor. 1699 my_wr=-999 1700 do io=1,nomegasf 1701 if (omegasf(io)>my_max_rest) then 1702 my_wr=io; EXIT 1703 end if 1704 end do 1705 if (my_wr==nomegasf+2) my_wr=nomegasf+1 1706 my_wl=-999 1707 do io=nomegasf,1,-1 1708 if (omegasf(io)< my_min_rest) then ! Check metals 1709 my_wl=io; EXIT 1710 end if 1711 end do 1712 1713 write(msg,'(a,2(1x,i0))')' my_wl and my_wr:',my_wl,my_wr 1714 call wrtout(std_out, msg) 1715 1716 if (my_wl==-999 .or. my_wr==-999) then 1717 write(msg,'(a,2i6)')' wrong value in my_wl and/or my_wr ',my_wl,my_wr 1718 ABI_ERROR(msg) 1719 end if 1720 1721 ! Calculate weights for Hilbert transform. 1722 call calc_kkweight(nomega,omega,nomegasf,omegasf,zcut,max_rest,kkweight) 1723 1724 end subroutine setup_spectral
m_chi0tk/symmetrize_afm_chi0tk [ Functions ]
[ Top ] [ m_chi0tk ] [ Functions ]
NAME
symmetrize_afm_chi0
FUNCTION
Reconstruct the (down, down) component of the irreducible polarizability starting from the (up,up) element in case of systems with AFM symmetries (i.e nspden==2 and nsppol=1). Return the trace (up,up)+(down,down) of the matrix as required by GW calculations.
INPUTS
Cryst<crystal_t>= Information on symmetries and unit cell. Gsph<gsphere_t>= The G-sphere used to descrive chi0. npwe=Number of G-vectors in chi0. nomega=number of frequencies. Ltg_q<littlegroup_t>=Structure with useful table describing the little group of the q-point.
SIDE EFFECTS
chi0(npwe,npwe,nomega)= In input the up-up component, in output the trace of chi0. The value of matrix elements that should be zero due to AFM symmetry properties are forced to be zero (see NOTES below). [chi0_lwing(npwe,nomega,3)] = Lower wings, symmetrized in output. [chi0_uwing(npwe,nomega,3)] = Upper wings, symmetrized in output. [chi0_head(3,3,nomega) ] = Head of chi0, symmetrized in output.
NOTES
In the case of magnetic group Shubnikov type III: For each set of paired FM-AFM symmetries, the down-down component of a generic response function in reciprocal space can be obtained according to: chi^{down,down}_{G1,G2}(q) = chi^{up up}_{G1,G2}(q) e^{iS(G1-G2).(tnonsFM - tnonsAFM)} where S is the rotational part common to the FM-AFM pair, tnonsFM and tnonsAFM are the fractional translations associated to the ferromagnetic and antiferromagnetic symmetry, respectively. Note that, if for a given G1-G2 pair, the phase e^{iS(G1-G2).(tnonsFM - tnonsAFM) depends on the FM-AFM symmetry pair, then the corresponding matrix element of chi0 must be zero. Actually this is manually enforced in the code because this property might not be perfectly satisfied due to round-off errors. In the case of magnetic group Shubnikov type III: Only the AFM symmetries that preserve the external q-point (with or without time-reversal) are used to get the (down, down) component using the fact that: chi^{down,down}_{G1,G2}(Sq) = chi^{up up}_{S^{-1}G1,S^{-1}G2}(q) e^{i(G2-G1).tnons_S } Actually we perform an average over subset of the little group of q with AFM character in order to reduce as much as possible errors due to round off errors. In brief we evaluate: 1/N_{Ltq} \sum_{S\in Ltg AFM} chi^{up up}_{S^{-1}G1,S^{-1}G2}(q) e^{i(G2-G1).tnons_S } where N_{Ltg} is the number of AFM operation in the little group (time reversal included)
TODO
It is possible to symmetrize chi0 without any the extra allocation for afm_mat. More CPU demanding but safer in case of a large chi0 matrix. One might loop over G1 and G2 shells ...
SOURCE
344 subroutine symmetrize_afm_chi0(Cryst,Gsph,Ltg_q,npwe,nomega,chi0,chi0_head,chi0_lwing,chi0_uwing) 345 346 !Arguments ------------------------------------ 347 !scalars 348 integer,intent(in) :: npwe,nomega 349 type(gsphere_t),intent(in) :: Gsph 350 type(crystal_t),intent(in) :: Cryst 351 type(littlegroup_t),intent(in) :: Ltg_q 352 !arrays 353 complex(gwpc),optional,intent(inout) :: chi0(npwe,npwe,nomega) 354 complex(dpc),optional,intent(inout) :: chi0_lwing(npwe,nomega,3) 355 complex(dpc),optional,intent(inout) :: chi0_uwing(npwe,nomega,3) 356 complex(dpc),optional,intent(inout) :: chi0_head(3,3,nomega) 357 358 !Local variables ------------------------------ 359 !scalars 360 integer :: io,ig1,ig2,isymf,isyma,isym,ipair,k0g,kg,npairs,nonzero 361 integer :: iSmg1,iSmg2,itim,shubnikov,ntest 362 complex(gwpc) :: phase,phase_old,sumchi,ctmp 363 logical :: found 364 !character(len=500) :: msg 365 !arrays 366 integer :: rotfm(3,3),rotafm(3,3),pairs2sym(2,Cryst%nsym/2) 367 real(dp) :: tfm(3),tafm(3) 368 complex(gwpc),allocatable :: afm_mat(:),chi0_afm(:,:) 369 370 !************************************************************************ 371 372 ABI_CHECK(ANY(Cryst%symafm==-1),'Not magnetic space group') 373 ! 374 ! ==== Find shubnikov type ==== 375 ! TODO This info should be stored in Cryst% 376 shubnikov=4; npairs=0 377 378 do isymf=1,Cryst%nsym 379 if (Cryst%symafm(isymf)==-1) CYCLE 380 rotfm = Cryst%symrec(:,:,isymf) 381 tfm = Cryst%tnons(:,isymf) 382 found = .FALSE. 383 384 do isyma=1,Cryst%nsym 385 if (Cryst%symafm(isyma)==1) CYCLE 386 rotafm = Cryst%symrec(:,:,isyma) 387 388 if (ALL(rotfm==rotafm)) then 389 found=.TRUE. 390 tafm = Cryst%tnons(:,isyma) 391 npairs=npairs+1 392 ABI_CHECK(npairs<=Cryst%nsym/2,'Wrong AFM group') 393 pairs2sym(1,npairs)=isymf 394 pairs2sym(2,npairs)=isyma 395 end if 396 end do !isyma 397 398 if (.not.found) then 399 shubnikov=3; EXIT !isymf 400 end if 401 end do !isymf 402 403 select case (shubnikov) 404 405 case (4) 406 call wrtout(std_out,' Found Magnetic group Shubnikov type IV') 407 ABI_CHECK(npairs==Cryst%nsym/2,'Wrong AFM space group') 408 409 ABI_MALLOC(afm_mat,(npwe*(npwe+1)/2)) 410 411 ! jmb 412 phase_old=zero 413 414 do ig2=1,npwe 415 k0g=ig2*(ig2-1)/2 416 do ig1=1,ig2 417 kg=k0g+ig1 418 nonzero=1 419 420 do ipair=1,Cryst%nsym/2 421 isymf = pairs2sym(1,ipair) 422 isyma = pairs2sym(2,ipair) 423 phase = ( Gsph%phmSGt(ig1,isymf)*conjg(Gsph%phmSGt(ig1,isyma)) ) * & 424 & ( Gsph%phmSGt(ig2,isymf)*conjg(Gsph%phmSGt(ig2,isyma)) ) 425 if (ipair>1 .and. (ABS(phase_old-phase) > tol6)) then 426 nonzero=0; EXIT 427 end if 428 phase_old=phase 429 end do !ipair 430 431 afm_mat(kg)=nonzero*(cone_gw + phase) 432 end do !ig1 433 end do !ig2 434 ! 435 ! ======================================================================= 436 ! ==== Symmetrize chi0 constructing chi0^{\up\up} + chi0{\down\down} ==== 437 ! ======================================================================= 438 ! 439 ! head^{\down\down} = head^{\up\up} 440 if (PRESENT(chi0_head)) chi0_head = two * chi0_head 441 ! 442 ! w^{\down\down}_{0 G'} = e^{-iSG'.(tFM-tAFM)} w^{\up\up}_{0 G'}. 443 ! w^{\down\down}_{G 0 } = e^{+iSG .(tFM-tAFM)} w^{\up\up}_{G 0 }. 444 if (PRESENT(chi0_uwing)) then 445 do io=1,nomega 446 do ig2=1,npwe 447 k0g=ig2*(ig2-1)/2 448 kg=k0g+1 449 chi0_uwing(ig2,io,:)=afm_mat(kg)*chi0_uwing(ig2,io,:) 450 end do 451 end do 452 end if 453 454 if (PRESENT(chi0_lwing)) then 455 do io=1,nomega 456 do ig1=1,npwe 457 k0g=ig1*(ig1-1)/2 458 kg=k0g+1 459 chi0_lwing(ig1,io,:)=conjg(afm_mat(kg))*chi0_lwing(ig1,io,:) 460 end do 461 end do 462 end if 463 464 if (PRESENT(chi0)) then 465 do io=1,nomega 466 ! Take care of diagonal. 467 do ig1=1,npwe 468 chi0(ig1,ig1,io)=two*chi0(ig1,ig1,io) 469 end do 470 471 ! Upper and lower triangle are treated differently: 472 ! We took advantage of the fact the afm_mat is hermitian to reduce memory. 473 do ig2=2,npwe 474 k0g=ig2*(ig2-1)/2 475 do ig1=1,ig2-1 476 kg=k0g+ig1 477 chi0(ig1,ig2,io)=afm_mat(kg)*chi0(ig1,ig2,io) 478 end do 479 end do 480 481 do ig1=2,npwe 482 k0g=ig1*(ig1-1)/2 483 do ig2=1,ig1-1 484 kg=k0g+ig2 485 chi0(ig1,ig2,io)=conjg(afm_mat(kg))*chi0(ig1,ig2,io) 486 end do 487 end do 488 end do !io 489 end if 490 491 ABI_FREE(afm_mat) 492 493 case (3) 494 call wrtout(std_out,' Found Magnetic group Shubnikov type III') 495 ABI_ERROR('Shubnikov type III not implemented') 496 497 ntest=0 498 do itim=1,ltg_q%timrev 499 do isym=1,ltg_q%nsym_sg 500 ! use only afm sym preserving q with and without time-reversal 501 if ( cryst%symafm(isym)==-1 .and. ltg_q%preserve(itim,isym)==1 ) ntest=ntest+1 502 end do 503 end do 504 505 if (ntest==0) then 506 ABI_WARNING("no symmetry can be used!") 507 end if 508 !RETURN 509 ABI_MALLOC(chi0_afm,(npwe,npwe)) 510 511 do io=1,nomega 512 513 do ig2=1,npwe 514 do ig1=1,npwe 515 sumchi=czero_gw 516 517 do itim=1,ltg_q%timrev 518 do isym=1,ltg_q%nsym_sg 519 ! use only afm sym preserving q with and without time-reversal 520 if ( cryst%symafm(isym)==-1 .and. ltg_q%preserve(itim,isym)==1 ) then 521 phase = Gsph%phmGt(ig1,isym)*conjg(Gsph%phmGt(ig2,isym)) 522 iSmg1=Gsph%rottbm1(ig1,itim,isym) 523 iSmg2=Gsph%rottbm1(ig2,itim,isym) 524 ctmp=chi0(iSmg1,iSmg2,io)*phase !; if (itim==2) ctmp=conjg(ctmp) !check this 525 sumchi=sumchi+ctmp !chi0(iSmg1,iSmg2,io)*phase 526 end if 527 end do ! isym 528 end do !itim 529 530 chi0_afm(ig1,ig2)=sumchi/Ltg_q%nsym_ltg !has to be changed in case of time-reversal 531 end do !ig1 532 end do !ig2 533 534 ! We want chi_{up,up} +chi_{dwn,dwn}. 535 chi0(:,:,io)=chi0(:,:,io)+chi0_afm(:,:) 536 end do !iomega 537 538 ABI_FREE(chi0_afm) 539 540 case default 541 ABI_BUG(sjoin('Wrong value for shubnikov= ', itoa(shubnikov))) 542 end select 543 544 end subroutine symmetrize_afm_chi0