TABLE OF CONTENTS
ABINIT/calc_sigc_cd [ Functions ]
NAME
calc_sigc_cd
FUNCTION
Calculate contributions to the self-energy operator with the contour deformation method.
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=Incremented with the number of poles whose contribution has not been taken into account due to limited frequency mesh used for W.
PARENTS
calc_sigc_me,m_screen
CHILDREN
spline,splint,xgemm,xgemv
SOURCE
45 #if defined HAVE_CONFIG_H 46 #include "config.h" 47 #endif 48 49 #include "abi_common.h" 50 51 52 subroutine calc_sigc_cd(npwc,npwx,nspinor,nomega,nomegae,nomegaer,nomegaei,rhotwgp,& 53 & omega,epsm1q,omegame0i,theta_mu_minus_e0i,ket,plasmafreq,npoles_missing,calc_poles,method) 54 55 use defs_basis 56 use m_profiling_abi 57 use m_splines 58 use m_xomp 59 60 use m_gwdefs, only : czero_gw, cone_gw, Kron15N, Kron15W, Gau7W, & 61 & Kron23N, Kron23W, Gau11W, Kron31N, Kron31W, Gau15W 62 use m_blas, only : xgemv, xgemm 63 64 !This section has been created automatically by the script Abilint (TD). 65 !Do not modify the following lines by hand. 66 #undef ABI_FUNC 67 #define ABI_FUNC 'calc_sigc_cd' 68 !End of the abilint section 69 70 implicit none 71 72 !Arguments ------------------------------------ 73 !scalars 74 integer,intent(in) :: nomega,nomegae,nomegaei,nomegaer,npwc,npwx,nspinor 75 integer,intent(inout) :: npoles_missing 76 real(dp),intent(in) :: theta_mu_minus_e0i,plasmafreq 77 !arrays 78 real(dp),intent(in) :: omegame0i(nomega) 79 complex(dpc),intent(in) :: omega(nomegae) 80 complex(gwpc),intent(in) :: epsm1q(npwc,npwc,nomegae) 81 complex(gwpc),intent(in) :: rhotwgp(npwx*nspinor) 82 complex(gwpc),intent(inout) :: ket(nspinor*npwc,nomega) 83 logical, intent(in), optional :: calc_poles(nomega) 84 integer, intent(in), optional :: method 85 86 !Local variables------------------------------- 87 !scalars 88 integer, parameter :: FABIEN=1,TRAPEZOID=2,NSPLINE=3 89 integer :: ii,ig,io,ios,ispinor,spadc,spadx,my_err,ierr,GK_LEVEL,INTMETHOD 90 integer :: i,j 91 real(dp) :: rt_imag,rt_real,local_one,local_zero 92 real(dp) :: intsign,temp1,temp2,temp3,temp4 93 real(dp) :: alph,inv_alph,beta,alphsq,betasq,inv_beta 94 real(dp) :: re_intG,re_intK,im_intG,im_intK,GKttab,tau,ttil 95 real(dp) :: ref,imf,r,s,r2,s2 96 complex(dpc) :: ct,domegaleft,domegaright 97 complex(gwpc) :: fact 98 !arrays 99 real(dp) :: omegame0i_tmp(nomega),tmp_x(2),tmp_y(2) 100 real(dp) :: left(nomega),right(nomega) 101 real(dp) :: tbeta(nomega),tinv_beta(nomega),tbetasq(nomega) 102 real(dp) :: atermr(nomega),aterml(nomega),logup(nomega),logdown(nomega) 103 real(dp) :: rtmp_r(nomegaer),rtmp_i(nomegaer) 104 real(dp) :: ftab(nomegaei+2),ftab2(nomegaei+2),xtab(nomegaei+2),y(3,nomegaei+2) 105 real(dp) :: work(nomegaei+2),work2(nomegaei+2),y2(3,nomegaei+2) 106 complex(dpc) :: omega_imag(nomegaei+1) 107 complex(gwpc) :: epsrho(npwc,nomegae),epsrho_imag(npwc,nomegaei+1) 108 complex(gwpc) :: tfone(npwc,nomegaei+1),tftwo(npwc,nomegaei+1) 109 complex(gwpc) :: weight(nomegaei+1,nomega) 110 complex(gwpc) :: weight2(nomegaei,nomega) 111 logical :: my_calc_poles(nomega) 112 real(dp), allocatable :: KronN(:),KronW(:),GaussW(:),fint(:),fint2(:) 113 !************************************************************************* 114 115 my_calc_poles=.TRUE.; my_err=0 116 117 ! Set integration method for imaginary axis 118 INTMETHOD = FABIEN 119 if (present(method)) then 120 if (method==1) INTMETHOD = FABIEN 121 if (method==2) INTMETHOD = TRAPEZOID 122 if (method>2) then 123 INTMETHOD = NSPLINE 124 if (method==3) then 125 GK_LEVEL = 15 126 ABI_ALLOCATE(KronN,(GK_LEVEL)) 127 ABI_ALLOCATE(KronW,(GK_LEVEL)) 128 ABI_ALLOCATE(GaussW,(GK_LEVEL-8)) 129 ABI_ALLOCATE(fint,(GK_LEVEL)) 130 ABI_ALLOCATE(fint2,(GK_LEVEL)) 131 KronN(:) = Kron15N(:); KronW(:) = Kron15W(:); GaussW(:) = Gau7W(:) 132 else if (method==4) then 133 GK_LEVEL = 23 134 ABI_ALLOCATE(KronN,(GK_LEVEL)) 135 ABI_ALLOCATE(KronW,(GK_LEVEL)) 136 ABI_ALLOCATE(GaussW,(GK_LEVEL-12)) 137 ABI_ALLOCATE(fint,(GK_LEVEL)) 138 ABI_ALLOCATE(fint2,(GK_LEVEL)) 139 KronN(:) = Kron23N(:); KronW(:) = Kron23W(:); GaussW(:) = Gau11W(:) 140 else if (method>4) then 141 GK_LEVEL = 31 142 ABI_ALLOCATE(KronN,(GK_LEVEL)) 143 ABI_ALLOCATE(KronW,(GK_LEVEL)) 144 ABI_ALLOCATE(GaussW,(GK_LEVEL-16)) 145 ABI_ALLOCATE(fint,(GK_LEVEL)) 146 ABI_ALLOCATE(fint2,(GK_LEVEL)) 147 KronN(:) = Kron31N(:); KronW(:) = Kron31W(:); GaussW(:) = Gau15W(:) 148 end if 149 end if 150 end if 151 152 ! Avoid divergences in $\omega - \omega_s$. 153 omegame0i_tmp(:)=omegame0i(:) 154 do ios=1,nomega 155 if (ABS(omegame0i_tmp(ios))<tol6) omegame0i_tmp(ios)=sign(tol6,omegame0i_tmp(ios)) 156 end do 157 158 do ispinor=1,nspinor 159 spadx=(ispinor-1)*npwx; spadc=(ispinor-1)*npwc 160 161 ! Calculate $ \sum_{Gp} (\epsilon^{-1}_{G Gp}(\omega)-\delta_{G Gp}) \rhotwgp(Gp) $ 162 !$omp parallel do 163 do io=1,nomegae 164 call XGEMV('N',npwc,npwc,cone_gw,epsm1q(:,:,io),npwc,rhotwgp(1+spadx:),1,czero_gw,epsrho(:,io),1) 165 end do 166 167 ! Integrand along the imaginary axis. 168 epsrho_imag(:,1)=epsrho(:,1) 169 epsrho_imag(:,2:nomegaei+1)=epsrho(:,nomegaer+1:nomegae) 170 171 ! Frequency mesh for integral along the imaginary axis. 172 omega_imag(1)=omega(1) 173 omega_imag(2:nomegaei+1)=omega(nomegaer+1:nomegae) 174 175 ! Original implementation -- saved here for reference during development 176 ! === Perform integration along the imaginary axis === 177 !do io=1,nomegaei+1 178 ! if (io==1) then 179 ! domegaleft = omega_imag(io) 180 ! domegaright =(omega_imag(io+1)-omega_imag(io ))*half 181 ! else if (io==nomegaei+1) then 182 ! domegaleft =(omega_imag(io )-omega_imag(io-1))*half 183 ! domegaright =(omega_imag(io )-omega_imag(io-1))*half 184 ! else 185 ! domegaleft =(omega_imag(io )-omega_imag(io-1))*half 186 ! domegaright =(omega_imag(io+1)-omega_imag(io ))*half 187 ! end if 188 ! do ios=1,nomega 189 ! omg2 = -AIMAG(omega_imag(io)+domegaright)/REAL(omegame0i_tmp(ios)) 190 ! omg1 = -AIMAG(omega_imag(io)-domegaleft )/REAL(omegame0i_tmp(ios)) 191 ! fact = ATAN(omg2)-ATAN(omg1) 192 ! ket(spadc+1:spadc+npwc,ios)=ket(spadc+1:spadc+npwc,ios)+epsrho_imag(:,io)*fact 193 ! end do 194 !end do !io 195 196 !ket(spadc+1:spadc+npwc,:)=ket(spadc+1:spadc+npwc,:)/pi 197 ! ---------------- end of original implementation ----------------------- 198 199 select case(INTMETHOD) 200 case(FABIEN) 201 ! Hopefully more effective implementation MS 04.08.2011 202 ! Perform integration along imaginary axis using BLAS 203 ! First calculate first and last weights 204 weight(1,:) = ATAN(-half*AIMAG(omega_imag(2))/REAL(omegame0i_tmp(:))) 205 domegaleft = (three*omega_imag(nomegaei+1)-omega_imag(nomegaei)) 206 domegaright = (omega_imag(nomegaei+1)+omega_imag(nomegaei)) 207 right(:) = -AIMAG(omega_imag(nomegaei+1)-omega_imag(nomegaei))*REAL(omegame0i_tmp(:)) 208 left(:) = quarter*AIMAG(domegaleft)*AIMAG(domegaright) & 209 & +REAL(omegame0i_tmp(:))*REAL(omegame0i_tmp(:)) 210 weight(nomegaei+1,:) = ATAN(right(:)/left(:)) 211 ! Calculate the rest of the weights 212 do io=2,nomegaei 213 domegaleft = (omega_imag(io )+omega_imag(io-1)) 214 domegaright = (omega_imag(io+1)+omega_imag(io )) 215 right(:) = -half*AIMAG(omega_imag(io+1)-omega_imag(io-1))*REAL(omegame0i_tmp(:)) 216 left(:) = REAL(omegame0i_tmp(:))*REAL(omegame0i_tmp(:)) & 217 & +quarter*AIMAG(domegaleft)*AIMAG(domegaright) 218 weight(io,:) = ATAN(right(:)/left(:)) 219 end do 220 221 ! Use BLAS call to perform matrix-matrix multiplication and accumulation 222 fact = CMPLX(piinv,zero) 223 224 call xgemm('N','N',npwc,nomega,nomegaei+1,fact,epsrho_imag,npwc,& 225 & weight,nomegaei+1,cone_gw,ket(spadc+1:spadc+npwc,:),npwc) 226 227 case(TRAPEZOID) 228 ! Trapezoidal rule Transform omega coordinates 229 alph = plasmafreq 230 alphsq = alph*alph 231 inv_alph = one/alph 232 233 xtab(1:nomegaei+1) = AIMAG(omega_imag(:))/(AIMAG(omega_imag(:)) + alph) 234 xtab(nomegaei+2) = one 235 236 ! Efficient trapezoidal rule with BLAS calls 237 tbeta(:) = REAL(omegame0i_tmp(:)) 238 tbetasq(:) = tbeta(:)*tbeta(:) 239 tinv_beta(:) = one/tbeta(:) 240 241 do io=1,nomegaei 242 atermr(:) = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(io+1)-tbetasq(:)) 243 aterml(:) = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(io )-tbetasq(:)) 244 right(:) = ATAN((atermr(:)-aterml(:))/(one+atermr(:)*aterml(:))) 245 logup(:) = ABS(((alphsq+tbetasq(:))*xtab(io+1)-two*tbetasq(:)) & 246 & *xtab(io+1)+tbetasq(:)) 247 logdown(:) = ABS(((alphsq+tbetasq(:))*xtab(io )-two*tbetasq(:)) & 248 & *xtab(io )+tbetasq(:)) 249 ! Trapezoid integration weights 250 weight(io,:) = CMPLX(-(half*alph*tbeta(:)*LOG(logup(:)/logdown(:)) + tbetasq(:) & 251 & *right(:))/(alphsq+tbetasq(:)),zero) 252 weight2(io,:) = CMPLX(-right(:),zero) 253 ! Linear interpolation coefficients for each section (sum over ig) 254 tfone(:,io) = (epsrho_imag(:,io+1)-epsrho_imag(:,io)) & 255 & /(xtab(io+1)-xtab(io)) 256 tftwo(:,io) = epsrho_imag(:,io) - tfone(:,io)*xtab(io) 257 end do 258 259 ! Calculate weights for asymptotic behaviour 260 atermr(:) = alph*tinv_beta(:) 261 aterml(:) = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(nomegaei+1)-tbetasq(:)) 262 logup(:) = alphsq*xtab(nomegaei+1)*xtab(nomegaei+1) 263 logdown(:) = ABS(((alphsq+tbetasq(:))*xtab(nomegaei+1)-two*tbetasq(:)) & 264 & *xtab(nomegaei+1)+tbetasq(:)) 265 right(:) = ATAN((atermr(:)-aterml(:))/(one+atermr(:)*aterml(:))) 266 weight (nomegaei+1,:) = CMPLX(-(half*(alphsq*tinv_beta(:)*LOG(logdown(:)/logup(:)) & 267 & - tbeta(:)*LOG(xtab(nomegaei+1)*xtab(nomegaei+1))) - alph*right(:)),zero) 268 tfone(:,nomegaei+1) = -(zero-epsrho_imag(:,nomegaei+1)*AIMAG(omega_imag(nomegaei+1))) & 269 /(one-xtab(nomegaei+1)) 270 271 ! Use BLAS call to perform matrix-matrix multiplication and accumulation 272 fact = CMPLX(piinv,zero) 273 274 call xgemm('N','N',npwc,nomega,nomegaei+1,fact,tfone,npwc,& 275 & weight ,nomegaei+1,cone_gw,ket(spadc+1:spadc+npwc,:),npwc) 276 call xgemm('N','N',npwc,nomega,nomegaei ,fact,tftwo,npwc,& 277 & weight2,nomegaei ,cone_gw,ket(spadc+1:spadc+npwc,:),npwc) 278 279 case(NSPLINE) 280 ! Natural spline followed by Gauss-Kronrod 281 ! Transform omega coordinates 282 alph = plasmafreq 283 alphsq = alph*alph 284 inv_alph = one/alph 285 286 xtab(1:nomegaei+1) = AIMAG(omega_imag(:))/(AIMAG(omega_imag(:)) + alph) 287 xtab(nomegaei+2) = one 288 289 ! Gauss-Kronrod integration of spline fit of f(t)/(1-t) in transformed space 290 ! *** OPENMP SECTION *** Added by MS 291 !!$OMP PARALLEL DO PRIVATE(ig,ftab,ftab2,s,s2,r,r2,y,y2,work,work2,beta,betasq,inv_beta, & 292 !!$OMP intsign,io,ii,i,j,re_intG,re_intK,im_intG,im_intK,temp1,temp2,temp3,temp4, & 293 !!$OMP ttil,tau,ref,fint,imf,fint2,GKttab) 294 do ig=1,npwc 295 ! Spline fit 296 ftab (1:nomegaei+1) = REAL(epsrho_imag(ig,1:nomegaei+1))/(one-xtab(1:nomegaei+1)) 297 ftab2(1:nomegaei+1) = AIMAG(epsrho_imag(ig,1:nomegaei+1))/(one-xtab(1:nomegaei+1)) 298 ftab (nomegaei+2) = zero; ftab2(nomegaei+2) = zero 299 ! Explicit calculation of spline coefficients 300 s = zero; s2 = zero 301 do i = 1, nomegaei+2-1 302 r = ( ftab (i+1) - ftab (i) ) / ( xtab(i+1) - xtab(i) ) 303 r2 = ( ftab2(i+1) - ftab2(i) ) / ( xtab(i+1) - xtab(i) ) 304 y (2,i) = r - s; y2(2,i) = r2 - s2 305 s = r; s2 = r2 306 end do 307 s = zero; s2 = zero 308 r = zero; r2 = zero 309 y(2,1) = zero; y2(2,1) = zero 310 y(2,nomegaei+2) = zero; y2(2,nomegaei+2) = zero 311 do i = 2, nomegaei+2-1 312 y (2,i) = y (2,i) + r * y (2,i-1) 313 y2(2,i) = y2(2,i) + r2 * y2(2,i-1) 314 work (i) = two * ( xtab(i-1) - xtab(i+1) ) - r * s 315 work2(i) = two * ( xtab(i-1) - xtab(i+1) ) - r2 * s2 316 s = xtab(i+1) - xtab(i) 317 s2 = s 318 r = s / work (i) 319 r2 = s2 / work2(i) 320 end do 321 do j = 2, nomegaei+2-1 322 i = nomegaei+2+1-j 323 y (2,i) = ( ( xtab(i+1) - xtab(i) ) * y (2,i+1) - y (2,i) ) / work (i) 324 y2(2,i) = ( ( xtab(i+1) - xtab(i) ) * y2(2,i+1) - y2(2,i) ) / work2(i) 325 end do 326 do i = 1, nomegaei+2-1 327 s = xtab(i+1) - xtab(i); s2 = s; 328 r = y(2,i+1) - y(2,i); r2 = y2(2,i+1) - y2(2,i); 329 y(3,i) = r / s; y2(3,i) = r2 / s2; 330 y(2,i) = three * y(2,i); y2(2,i) = three * y2(2,i); 331 y (1,i) = ( ftab (i+1) - ftab (i) ) / s - ( y (2,i) + r ) * s 332 y2(1,i) = ( ftab2(i+1) - ftab2(i) ) / s2 - ( y2(2,i) + r2 ) * s2 333 end do 334 ! End of spline interpolation 335 do ios=1,nomega 336 beta = REAL(omegame0i_tmp(ios)) 337 betasq = beta*beta 338 inv_beta = one/beta 339 intsign = sign(half*piinv,beta) 340 beta = ABS(beta) 341 io = 1; re_intG = zero; re_intK = zero; im_intG = zero; im_intK = zero 342 do ii=1,GK_LEVEL 343 do 344 GKttab = two*alph*xtab(io+1)/(beta-(beta-alph)*xtab(io+1))-one 345 if (GKttab > KronN(ii)) EXIT 346 io = io + 1 347 end do 348 temp1 = half*(KronN(ii)+one) 349 temp2 = temp1 - half 350 temp3 = temp2*temp2 351 temp4 = half/(temp3 + quarter) 352 ttil = beta*temp1/(alph-(alph-beta)*temp1) 353 tau = ttil - xtab(io) 354 ref = ftab (io) + tau*(y (1,io)+tau*(y (2,io)+tau*y (3,io))) 355 fint (ii) = -ref*(one-ttil)*temp4 356 imf = ftab2(io) + tau*(y2(1,io)+tau*(y2(2,io)+tau*y2(3,io))) 357 fint2(ii) = -imf*(one-ttil)*temp4 358 re_intK = KronW(ii)*fint (ii) 359 im_intK = KronW(ii)*fint2(ii) 360 ket(spadc+ig,ios) = ket(spadc+ig,ios)+intsign*CMPLX(re_intK,im_intK) 361 end do ! ii 362 end do !ios 363 end do !ig 364 !!$OMP END PARALLEL DO 365 366 end select !intmethod 367 368 local_one = one 369 local_zero = zero 370 371 ! ============================================ 372 ! ==== Add contribution coming from poles ==== 373 ! ============================================ 374 ! First see if the contribution has been checked before the routine is entered 375 if (present(calc_poles)) then 376 my_calc_poles = calc_poles 377 else ! Otherwise check locally if there is a contribution 378 do ios=1,nomega 379 if (omegame0i_tmp(ios)>tol12) then 380 if ((local_one-theta_mu_minus_e0i)<tol12) my_calc_poles(ios) = .FALSE. 381 end if 382 if (omegame0i_tmp(ios)<-tol12) then 383 if (theta_mu_minus_e0i<tol12) my_calc_poles(ios) = .FALSE. 384 end if 385 end do !ios 386 end if 387 388 if (ANY(my_calc_poles(:))) then ! Make sure we only enter if necessary 389 ! *** OPENMP SECTION *** Added by MS 390 !!OMP !write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',xomp_get_num_threads() 391 !$OMP PARALLEL SHARED(npwc,nomega,nomegaer,theta_mu_minus_e0i,spadc,local_one,local_zero, & 392 !$OMP omega,epsrho,omegame0i_tmp,ket,my_calc_poles) & 393 !$OMP PRIVATE(ig,ios,rtmp_r,rtmp_i,tmp_x,tmp_y,rt_real,rt_imag,ct,ierr) REDUCTION(+:my_err) 394 !!OMP $ write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',xomp_get_num_threads() 395 !$OMP DO 396 do ig=1,npwc 397 ! Prepare the spline interpolation by filling at once the arrays rtmp_r, rtmp_i 398 call spline(DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),nomegaer,local_zero,local_zero,rtmp_r) 399 call spline(DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),nomegaer,local_zero,local_zero,rtmp_i) 400 ! call spline_complex( DBLE(omega(1:nomegaer)), epsrho(ig,1:nomegaer), nomegaer, zero, zero, rtmp ) 401 402 do ios=1,nomega 403 if (.NOT.my_calc_poles(ios)) CYCLE 404 405 ! Interpolate real and imaginary part of epsrho at |omegame0i_tmp|. 406 tmp_x(1) = ABS(omegame0i_tmp(ios)) 407 call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),rtmp_r,1,tmp_x,tmp_y,ierr=ierr) 408 if (ig==1.and.ispinor==1) my_err = my_err + ierr 409 rt_real = tmp_y(1) 410 411 tmp_x(1) = ABS(omegame0i_tmp(ios)) 412 call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),rtmp_i,1,tmp_x,tmp_y) 413 rt_imag = tmp_y(1) 414 !!call splint_complex(nomegaer,DBLE(omega(1:nomegaer)),epsrho(ig,1:nomegaer),rtmp,1,tmp_x,yfit) 415 416 ct=DCMPLX(rt_real,rt_imag) 417 418 if (omegame0i_tmp(ios)>tol12) then 419 ket(spadc+ig,ios)=ket(spadc+ig,ios)+ct*(local_one-theta_mu_minus_e0i) 420 end if 421 if (omegame0i_tmp(ios)<-tol12) then 422 ket(spadc+ig,ios)=ket(spadc+ig,ios)-ct*theta_mu_minus_e0i 423 end if 424 425 end do !ios 426 end do !ig 427 !$OMP END DO 428 !$OMP END PARALLEL 429 end if ! ANY(my_calc_poles) 430 end do !ispinor 431 432 npoles_missing = npoles_missing + my_err 433 434 if (INTMETHOD>2) then 435 ABI_DEALLOCATE(KronN) 436 ABI_DEALLOCATE(KronW) 437 ABI_DEALLOCATE(GaussW) 438 ABI_DEALLOCATE(fint) 439 ABI_DEALLOCATE(fint2) 440 end if 441 442 end subroutine calc_sigc_cd