TABLE OF CONTENTS
ABINIT/calc_sigc_pole_cd [ Functions ]
NAME
calc_sigc_pole_cd
FUNCTION
Calculate contributions to the self-energy operator with the contour deformation method, using a pole-fit screening.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (FB, GMR, VO, LR, RWG, MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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=Incremeented with the number of poles whose contribution has not been taken into account due to limited frequency mesh used for W.
PARENTS
CHILDREN
cgemm,re_and_im_screening_with_phase,spline,splint,xgemv,zgemm
SOURCE
45 #if defined HAVE_CONFIG_H 46 #include "config.h" 47 #endif 48 49 #include "abi_common.h" 50 51 52 subroutine calc_sigc_pole_cd(npwc,npwx,nspinor,ncoeff,nomega,nomegae,nomegaer,nomegaei,rhotwgp,& 53 & omega,epsm1q,omegame0i,theta_mu_minus_e0i,ket,npoles_missing,calc_poles) 54 55 use m_profiling_abi 56 57 use defs_basis 58 use m_splines 59 60 use m_gwdefs, only : czero_gw, cone_gw 61 use m_blas, only : xgemv 62 use m_model_screening, only : re_and_im_screening_with_phase 63 !$ use omp_lib, only: omp_get_thread_num, omp_get_num_threads 64 65 !This section has been created automatically by the script Abilint (TD). 66 !Do not modify the following lines by hand. 67 #undef ABI_FUNC 68 #define ABI_FUNC 'calc_sigc_pole_cd' 69 !End of the abilint section 70 71 implicit none 72 73 !Arguments ------------------------------------ 74 !scalars 75 integer,intent(in) :: nomega,ncoeff,nomegae,nomegaei,nomegaer,npwc,npwx,nspinor 76 integer,intent(inout) :: npoles_missing 77 real(dp),intent(in) :: theta_mu_minus_e0i 78 !arrays 79 real(dp),intent(in) :: omegame0i(nomega) 80 complex(dpc),intent(in) :: omega(nomegae) 81 real(gwp) :: epsm1q(npwc,npwc,ncoeff) 82 complex(gwpc),intent(in) :: rhotwgp(npwx*nspinor) 83 complex(gwpc),intent(inout) :: ket(nspinor*npwc,nomega) 84 logical, intent(in), optional :: calc_poles(nomega) 85 86 !Local variables------------------------------- 87 !scalars 88 integer :: ig,io,ios,ispinor,spadc,spadx,my_err,ig1,ig2 89 real(dp) :: rt_imag,rt_real,local_one,local_zero 90 ! real(dp) :: intresult,intsign 91 complex(dpc) :: ct,domegaleft,domegaright 92 #if defined HAVE_GW_DPC 93 complex(dpc) :: fact 94 #else 95 complex(spc) :: fact 96 #endif 97 !arrays 98 real(dp) :: omegame0i_tmp(nomega),tmp_x(2),tmp_y(2) 99 real(dp) :: rtmp_r(nomegaer),rtmp_i(nomegaer) 100 ! real(dp) :: ftab(nomegaei+2),xtab(nomegaei+2),y(3,nomegaei+2) 101 ! real(dp) :: work(nomegaei+2),e(nomegaei+2) 102 complex(dpc) :: omega_imag(nomegaei+1) 103 complex(gwpc) :: epsrho(npwc,nomegae),epsrho_imag(npwc,nomegaei+1) 104 complex(gwpc) :: weight(nomegaei+1,nomega) 105 complex(gwpc) :: work(npwc,npwc,nomegae) 106 logical :: my_calc_poles(nomega) 107 !************************************************************************* 108 109 my_err=0 110 my_calc_poles=.TRUE. 111 112 ! Avoid divergences in $\omega - \omega_s$. 113 omegame0i_tmp(:)=omegame0i(:) 114 do ios=1,nomega 115 if (ABS(omegame0i_tmp(ios))<tol6) omegame0i_tmp(ios)=tol6 116 end do 117 118 do ispinor=1,nspinor 119 spadx=(ispinor-1)*npwx 120 spadc=(ispinor-1)*npwc 121 ! 122 ! Calculate $ \sum_{Gp} (\epsilon^{-1}_{G Gp}(\omega)-\delta_{G Gp}) \rhotwgp(Gp) $ 123 ! First calculate the values from the poles 124 do ig2=1,npwc 125 do ig1=1,npwc 126 call re_and_im_screening_with_phase(omega,work(ig1,ig2,:),nomegae, & 127 & epsm1q(ig1,ig2,:),ncoeff) 128 end do 129 end do 130 131 do io=1,nomegae 132 call XGEMV('N',npwc,npwc,cone_gw,work(:,:,io),npwc,rhotwgp(1+spadx:),1,czero_gw,epsrho(:,io),1) 133 end do 134 135 !write(std_out,*) ch10,' Epsrho: ',epsrho 136 137 ! Integrand along the imaginary axis. 138 epsrho_imag(:,1)=epsrho(:,1) 139 epsrho_imag(:,2:nomegaei+1)=epsrho(:,nomegaer+1:nomegae) 140 141 ! Frequency mesh for integral along the imaginary axis. 142 omega_imag(1)=omega(1) 143 omega_imag(2:nomegaei+1)=omega(nomegaer+1:nomegae) 144 145 ! Original implementation -- saved here for reference during development 146 ! === Perform integration along the imaginary axis === 147 !do io=1,nomegaei+1 148 ! 149 ! if (io==1) then 150 ! domegaleft = omega_imag(io) 151 ! domegaright =(omega_imag(io+1)-omega_imag(io ))*half 152 ! else if (io==nomegaei+1) then 153 ! domegaleft =(omega_imag(io )-omega_imag(io-1))*half 154 ! domegaright =(omega_imag(io )-omega_imag(io-1))*half 155 ! else 156 ! domegaleft =(omega_imag(io )-omega_imag(io-1))*half 157 ! domegaright =(omega_imag(io+1)-omega_imag(io ))*half 158 ! end if 159 ! do ios=1,nomega 160 ! omg2 = -AIMAG(omega_imag(io)+domegaright)/REAL(omegame0i_tmp(ios)) 161 ! omg1 = -AIMAG(omega_imag(io)-domegaleft )/REAL(omegame0i_tmp(ios)) 162 ! fact = ATAN(omg2)-ATAN(omg1) 163 ! ket(spadc+1:spadc+npwc,ios)=ket(spadc+1:spadc+npwc,ios)+epsrho_imag(:,io)*fact 164 ! end do 165 !end do !io 166 167 !ket(spadc+1:spadc+npwc,:)=ket(spadc+1:spadc+npwc,:)/pi 168 ! 169 ! ---------------- end of original implementation ----------------------- 170 171 ! Hopefully more effective implementation MS 04.08.2011 172 ! === Perform integration along imaginary axis using BLAS === 173 ! First calculate first and last weights 174 domegaleft = omega_imag(1) 175 domegaright = (omega_imag(2)-omega_imag(1))*half 176 weight(1,:) = ATAN(-AIMAG(omega_imag(1)+domegaright)/REAL(omegame0i_tmp(:)))& 177 & -ATAN(-AIMAG(omega_imag(1)-domegaleft)/REAL(omegame0i_tmp(:))) 178 domegaleft = (omega_imag(nomegaei+1)-omega_imag(nomegaei))*half 179 domegaright = (omega_imag(nomegaei+1)-omega_imag(nomegaei))*half 180 weight(nomegaei+1,:) = ATAN(-AIMAG(omega_imag(nomegaei+1)+domegaright)/REAL(omegame0i_tmp(:)))& 181 & -ATAN(-AIMAG(omega_imag(nomegaei+1)-domegaleft)/REAL(omegame0i_tmp(:))) 182 ! Calculate the rest of the weights 183 do io=2,nomegaei 184 domegaleft = (omega_imag(io )-omega_imag(io-1))*half 185 domegaright = (omega_imag(io+1)-omega_imag(io ))*half 186 weight(io,:) = ATAN(-AIMAG(omega_imag(io)+domegaright)/REAL(omegame0i_tmp(:)))& 187 & -ATAN(-AIMAG(omega_imag(io)-domegaleft)/REAL(omegame0i_tmp(:))) 188 end do 189 190 ! Use BLAS call to perform matrix-matrix multiplication and accumulation 191 fact = CMPLX(piinv,zero) 192 #if defined HAVE_GW_DPC 193 call zgemm('N','N',npwc,nomega,nomegaei+1,fact,epsrho_imag,npwc,& 194 & weight,nomegaei+1,cone,ket(spadc+1:spadc+npwc,:),npwc) 195 #else 196 call cgemm('N','N',npwc,nomega,nomegaei+1,fact,epsrho_imag,npwc,& 197 & weight,nomegaei+1,CMPLX(1.0_sp,0.0_sp),ket(spadc+1:spadc+npwc,:),npwc) 198 #endif 199 200 !#else 201 ! Use spline integration 202 !do ios=1,nomega 203 ! if (ABS(omegame0i(ios))<tol12) then 204 ! intsign = sign(one,-REAL(omegame0i(ios))) 205 ! ket(spadc+1:spadc+npwc,ios) = ket(spadc+1:spadc+npwc,ios) + half*intsign*epsrho_imag(:,1) 206 ! else 207 ! xtab(1:nomegaei+1) = ABS(two*piinv*ATAN(-AIMAG(omega_imag(:))/REAL(omegame0i(ios)))) 208 ! intsign = sign(one,-REAL(omegame0i(ios))) 209 ! xtab(nomegaei+2) = one 210 ! do ig=1,npwc 211 ! ftab(1:nomegaei+1) = epsrho_imag(ig,:) 212 ! ftab(nomegaei+2) = zero 213 ! call cspint(ftab,xtab,nomegaei+2,zero,one,y,e,work,intresult) 214 ! ket(spadc+ig,ios) = ket(spadc+ig,ios) + half*intsign*intresult 215 ! end do 216 ! end if 217 !end do 218 !#endif 219 220 local_one = one 221 local_zero = zero 222 223 ! ============================================ 224 ! ==== Add contribution coming from poles ==== 225 ! ============================================ 226 ! First see if the contribution has been checked before the routine is entered 227 if (present(calc_poles)) then 228 my_calc_poles = calc_poles 229 else ! Otherwise check locally if there is a contribution 230 do ios=1,nomega 231 if (omegame0i_tmp(ios)>tol12) then 232 if ((local_one-theta_mu_minus_e0i)<tol12) my_calc_poles(ios) = .FALSE. 233 end if 234 if (omegame0i_tmp(ios)<-tol12) then 235 if (theta_mu_minus_e0i<tol12) my_calc_poles(ios) = .FALSE. 236 end if 237 end do !ios 238 end if 239 240 if (ANY(my_calc_poles(:))) then ! Make sure we only enter if necessary 241 242 !!OMP *** OPENMP SECTION *** Added by MS 243 !!OMP!$ write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',omp_get_num_threads() 244 !$OMP PARALLEL SHARED(npwc,nomega,nomegaer,theta_mu_minus_e0i,spadc,local_one,local_zero, & 245 !$OMP omega,epsrho,omegame0i_tmp,ket,my_calc_poles) & 246 !$OMP PRIVATE(ig,ios,rtmp_r,rtmp_i,tmp_x,tmp_y,rt_real,rt_imag,ct) 247 !!OMP!$ write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',omp_get_num_threads() 248 !$OMP DO 249 do ig=1,npwc 250 ! 251 ! * Prepare the spline interpolation by filling at once the arrays rtmp_r, rtmp_i 252 call spline(DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),nomegaer,local_zero,local_zero,rtmp_r) 253 call spline(DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),nomegaer,local_zero,local_zero,rtmp_i) 254 255 !! call spline_complex( DBLE(omega(1:nomegaer)), epsrho(ig,1:nomegaer), nomegaer, zero, zero, rtmp ) 256 257 do ios=1,nomega 258 259 if (.NOT.my_calc_poles(ios)) CYCLE 260 261 ! * Interpolate real and imaginary part of epsrho at |omegame0i_tmp|. 262 tmp_x(1) = ABS(omegame0i_tmp(ios)) 263 call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),rtmp_r,1,tmp_x,tmp_y,ierr=my_err) 264 if (ig==1.and.ispinor==1) npoles_missing = npoles_missing + my_err 265 rt_real = tmp_y(1) 266 267 tmp_x(1) = ABS(omegame0i_tmp(ios)) 268 call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),rtmp_i,1,tmp_x,tmp_y) 269 rt_imag = tmp_y(1) 270 271 !!call splint_complex(nomegaer,DBLE(omega(1:nomegaer)),epsrho(ig,1:nomegaer),rtmp,1,tmp_x,yfit) 272 273 ct=DCMPLX(rt_real,rt_imag) 274 275 if (omegame0i_tmp(ios)>tol12) then 276 ket(spadc+ig,ios)=ket(spadc+ig,ios)+ct*(local_one-theta_mu_minus_e0i) 277 end if 278 if (omegame0i_tmp(ios)<-tol12) then 279 ket(spadc+ig,ios)=ket(spadc+ig,ios)-ct*theta_mu_minus_e0i 280 end if 281 282 end do !ios 283 end do !ig 284 !$OMP END DO 285 !$OMP END PARALLEL 286 287 end if ! ANY(my_calc_poles) 288 289 end do !ispinor 290 291 end subroutine calc_sigc_pole_cd