TABLE OF CONTENTS
ABINIT/dieltcel [ Functions ]
NAME
dieltcel
FUNCTION
Compute either test charge or electronic dielectric matrices from susceptibility matrix Diagonalize it, then invert it.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, LSI) 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
gmet(3,3)=reciprocal space metric tensor in bohr**-2. kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix. kxc(nfft,nkxc)=exchange-correlation kernel, needed if the electronic dielectric matrix is computed nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkxc=second dimension of the array kxc, see rhotoxc.f for a description npwdiel=size of the dielinv and susmat arrays. nspden=number of spin-density components occopt=option for occupancies option=1 for Test Charge dielectric matrix, 2 for electronic dielectric matrix prtvol=control print volume and debugging output susmat(2,npwdiel,nspden,npwdiel,nspden)= the susceptibility (or density-density response) matrix in reciprocal space
OUTPUT
dielinv(2,npwdiel,nspden,npwdiel,nspden)=inverse of the (non-hermitian)
TC dielectric matrix in reciprocal space.
NOTES
Output (not cleaned) !!! Spin behaviour is not obvious !!! Will not work in the spin-polarized, metallic case.
PARENTS
prcref,prcref_PMA
CHILDREN
destroy_mpi_enreg,fourdp,init_distribfft_seq,initmpi_seq,timab,wrtout zhpev
SOURCE
52 #if defined HAVE_CONFIG_H 53 #include "config.h" 54 #endif 55 56 #include "abi_common.h" 57 58 subroutine dieltcel(dielinv,gmet,kg_diel,kxc,& 59 & nfft,ngfft,nkxc,npwdiel,nspden,occopt,option,paral_kgb,prtvol,susmat) 60 61 use defs_basis 62 use defs_abitypes 63 use m_errors 64 use m_profiling_abi 65 66 use m_mpinfo, only : destroy_mpi_enreg 67 use m_distribfft, only : init_distribfft_seq 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'dieltcel' 73 use interfaces_14_hidewrite 74 use interfaces_18_timing 75 use interfaces_51_manage_mpi 76 use interfaces_53_ffts 77 !End of the abilint section 78 79 implicit none 80 81 !Arguments ------------------------------------ 82 !scalars 83 integer,intent(in) :: nfft,nkxc,npwdiel,nspden,occopt,option,paral_kgb 84 integer,intent(in) :: prtvol 85 !arrays 86 integer,intent(in) :: kg_diel(3,npwdiel),ngfft(18) 87 real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc) 88 real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 89 real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden) 90 91 !Local variables------------------------------- 92 !scalars 93 integer :: i1,i2,i3,ieig,ier,ifft,ii,index,ipw0,ipw1,ipw2,ispden,j1 94 integer :: j2,j3,jj,k1,k2,k3,n1,n2,n3 95 real(dp) :: ai,ai2,ar,ar2,eiginv,gred1,gred2,gred3,gsquar,si 96 real(dp) :: sr,tpisq 97 character(len=500) :: message 98 type(MPI_type) :: mpi_enreg_seq 99 !arrays 100 real(dp) :: tsec(2) 101 real(dp),allocatable :: eig_msusinvsqr(:),eig_msussqr(:) 102 real(dp),allocatable :: eig_sus(:),eig_sym(:),invsqrsus(:,:,:,:,:) 103 real(dp),allocatable :: khxc(:,:,:,:,:),kxcg(:,:),sqrsus(:,:,:,:,:),sush(:) 104 real(dp),allocatable :: susvec(:,:,:),symdielmat(:,:,:,:,:),symh(:) 105 real(dp),allocatable :: symvec(:,:,:,:,:),wkxc(:),work(:,:,:,:,:) 106 real(dp),allocatable :: work2(:,:,:,:,:),zhpev1(:,:),zhpev2(:) 107 !no_abirules 108 !integer :: ipw3 109 !real(dp) :: elementi,elementr 110 !DEBUG 111 !Used to moderate divergence effect near rho=0 of the Kxc 112 !this limit value is truly empirical (exprmt on small Sr cell). 113 !real(dp) :: kxc_min=-200.0 114 !ENDDEBUG 115 116 ! ************************************************************************* 117 118 call timab(96,1,tsec) 119 120 !tpisq is (2 Pi) **2: 121 tpisq=(two_pi)**2 122 123 if(nspden/=1 .and. (occopt>=3 .and. occopt<=8) )then 124 write(message, '(a,a,a)' )& 125 & ' In the present version of the code, one cannot produce',ch10,& 126 & ' the dielectric matrix in the metallic, spin-polarized case.' 127 MSG_BUG(message) 128 end if 129 130 if(nspden==4)then 131 write(message,'(a,a,a)')& 132 & ' In the present version of the code, one cannot produce',ch10,& 133 & ' the dielectric matrix in the non-collinear spin-polarized case.' 134 MSG_ERROR(message) 135 end if 136 137 138 !-Diagonalize the susceptibility matrix 139 140 ABI_ALLOCATE(sush,(npwdiel*(npwdiel+1))) 141 ABI_ALLOCATE(susvec,(2,npwdiel,npwdiel)) 142 ABI_ALLOCATE(eig_msusinvsqr,(npwdiel)) 143 ABI_ALLOCATE(eig_msussqr,(npwdiel)) 144 ABI_ALLOCATE(eig_sus,(npwdiel)) 145 ABI_ALLOCATE(zhpev1,(2,2*npwdiel-1)) 146 ABI_ALLOCATE(zhpev2,(3*npwdiel-2)) 147 ABI_ALLOCATE(work,(2,npwdiel,nspden,npwdiel,nspden)) 148 ABI_ALLOCATE(work2,(2,npwdiel,nspden,npwdiel,nspden)) 149 ABI_ALLOCATE(sqrsus,(2,npwdiel,nspden,npwdiel,nspden)) 150 ABI_ALLOCATE(invsqrsus,(2,npwdiel,nspden,npwdiel,nspden)) 151 152 !At some time, should take care of different spin channels 153 do ispden=1,nspden 154 155 if(nspden/=1)then 156 MSG_ERROR('dieltcel : stop, nspden/=1') 157 end if 158 159 ! Store the susceptibility matrix in proper mode before calling zhpev 160 index=1 161 do ii=1,npwdiel 162 do jj=1,ii 163 sush(index )=susmat(1,jj,1,ii,1) 164 sush(index+1)=susmat(2,jj,1,ii,1) 165 index=index+2 166 end do 167 end do 168 169 ier=0 170 call ZHPEV ('V','U',npwdiel,sush,eig_sus,susvec,npwdiel,zhpev1,& 171 & zhpev2,ier) 172 173 ! DEBUG 174 ! write(std_out,*)' dieltcel : print eigenvalues of the susceptibility matrix' 175 ! do ii=1,npwdiel 176 ! write(std_out,'(i5,es16.6)' )ii,eig_sus(ii) 177 ! end do 178 ! ENDDEBUG 179 180 do ii=1,npwdiel 181 if(-eig_sus(ii)>1.d-12)then 182 eig_msussqr(ii)=sqrt(-eig_sus(ii)) 183 eig_msusinvsqr(ii)=1._dp/eig_msussqr(ii) 184 else if(-eig_sus(ii)< -1.d-12)then 185 message = "Found positive eigenvalue of susceptibility matrix." 186 MSG_BUG(message) 187 else 188 ! Set the eigenvalue corresponding to a constant potential change to 1, 189 ! while it will be set to zero in Khx. 190 eig_msussqr(ii)=1._dp 191 eig_msusinvsqr(ii)=1._dp 192 end if 193 end do 194 195 ! Compute square root of minus susceptibility matrix 196 ! and inverse square root of minus susceptibility matrix 197 do ii=1,npwdiel 198 work(:,:,1,ii,1)=susvec(:,:,ii)*eig_msussqr(ii) 199 work2(:,:,1,ii,1)=susvec(:,:,ii)*eig_msusinvsqr(ii) 200 end do 201 do ipw2=1,npwdiel 202 do ipw1=ipw2,npwdiel 203 ar=0._dp ; ai=0._dp ; ar2=0._dp ; ai2=0._dp 204 do ii=1,npwdiel 205 sr=susvec(1,ipw2,ii) ; si=susvec(2,ipw2,ii) 206 ar =ar +work(1,ipw1,1,ii,1)*sr +work(2,ipw1,1,ii,1)*si 207 ai =ai +work(2,ipw1,1,ii,1)*sr -work(1,ipw1,1,ii,1)*si 208 ar2=ar2 +work2(1,ipw1,1,ii,1)*sr +work2(2,ipw1,1,ii,1)*si 209 ai2=ai2 +work2(2,ipw1,1,ii,1)*sr -work2(1,ipw1,1,ii,1)*si 210 end do 211 sqrsus(1,ipw1,1,ipw2,1)=ar 212 sqrsus(2,ipw1,1,ipw2,1)=ai 213 invsqrsus(1,ipw1,1,ipw2,1)=ar2 214 invsqrsus(2,ipw1,1,ipw2,1)=ai2 215 if(ipw1/=ipw2)then 216 sqrsus(1,ipw2,1,ipw1,1)=ar 217 sqrsus(2,ipw2,1,ipw1,1)=-ai 218 invsqrsus(1,ipw2,1,ipw1,1)=ar2 219 invsqrsus(2,ipw2,1,ipw1,1)=-ai2 220 end if 221 end do 222 end do 223 224 ! DEBUG 225 ! Checks whether sqrsus and invsqrsus are inverse of each other. 226 ! do ipw1=1,npwdiel 227 ! do ipw2=1,npwdiel 228 ! elementr=0.0_dp 229 ! elementi=0.0_dp 230 ! do ipw3=1,npwdiel 231 ! elementr=elementr+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1)& 232 ! & -sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1) 233 ! elementi=elementi+sqrsus(1,ipw1,1,ipw3,1)*invsqrsus(2,ipw3,1,ipw2,1)& 234 ! & +sqrsus(2,ipw1,1,ipw3,1)*invsqrsus(1,ipw3,1,ipw2,1) 235 ! end do 236 ! if(elementr**2+elementi**2 > 1.0d-12)then 237 ! if( ipw1 /= ipw2 .or. & 238 ! & ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then 239 ! write(std_out,*)' dieltcel : sqrsus and invsqrsus are not (pseudo)',& 240 ! & 'inverse of each other' 241 ! write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2 242 ! write(std_out,*)' elementr,elementi=',elementr,elementi 243 ! stop 244 ! end if 245 ! end if 246 ! end do 247 ! end do 248 ! ENDDEBUG 249 250 ! End loop over spins 251 end do 252 253 ABI_DEALLOCATE(eig_msusinvsqr) 254 ABI_DEALLOCATE(eig_msussqr) 255 ABI_DEALLOCATE(eig_sus) 256 ABI_DEALLOCATE(sush) 257 ABI_DEALLOCATE(susvec) 258 259 !-Compute the Hxc kernel 260 261 ABI_ALLOCATE(khxc,(2,npwdiel,nspden,npwdiel,nspden)) 262 ABI_ALLOCATE(symdielmat,(2,npwdiel,nspden,npwdiel,nspden)) 263 264 khxc(:,:,:,:,:)=0.0_dp 265 266 !Compute Hartree kernel 267 do ipw1=1,npwdiel 268 gred1=dble(kg_diel(1,ipw1)) 269 gred2=dble(kg_diel(2,ipw1)) 270 gred3=dble(kg_diel(3,ipw1)) 271 gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 & 272 & +2.0_dp*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 + & 273 & gmet(2,3)*gred2*gred3) ) 274 ! Distinguish G=0 from other elements 275 if(gsquar>1.0d-12)then 276 khxc(1,ipw1,1,ipw1,1)= 4.0_dp*pi/gsquar 277 else 278 ! G=0 279 ipw0=ipw1 280 end if 281 end do 282 283 !Eventually add the xc part 284 if(option>=2)then 285 286 ABI_ALLOCATE(wkxc,(nfft)) 287 ABI_ALLOCATE(kxcg,(2,nfft)) 288 wkxc(:)=kxc(:,1) 289 ! DEBUG 290 ! Used to moderate divergenc effect near rho=0 of the Kxc (see above). 291 ! wkxc(:)=merge(kxc(:,1), kxc_min, kxc(:,1) > kxc_min) 292 ! ENDDEBUG 293 call initmpi_seq(mpi_enreg_seq) 294 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all') 295 call fourdp(1,kxcg,wkxc,-1,mpi_enreg_seq,nfft,ngfft,paral_kgb,0) ! trsfrm R to G 296 call destroy_mpi_enreg(mpi_enreg_seq) 297 298 ! Compute difference in G vectors 299 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 300 do ipw2=1,npwdiel 301 if(ipw2/=ipw0)then 302 303 j1=kg_diel(1,ipw2) ; j2=kg_diel(2,ipw2) ; j3=kg_diel(3,ipw2) 304 ! Fills diagonal 305 khxc(1,ipw2,1,ipw2,1)=khxc(1,ipw2,1,ipw2,1)+kxcg(1,1) 306 khxc(2,ipw2,1,ipw2,1)=khxc(2,ipw2,1,ipw2,1)+kxcg(2,1) 307 308 if(ipw2/=npwdiel)then 309 ! Fills off-diagonal part of the matrix, except G=0 310 do ipw1=ipw2+1,npwdiel 311 if(ipw1/=ipw0)then 312 i1=kg_diel(1,ipw1) ; i2=kg_diel(2,ipw1) ; i3=kg_diel(3,ipw1) 313 ! Use of two mod calls handles both i1-j1>=ndiel1 AND i1-j1<0 314 k1=mod(n1+mod(i1-j1,n1),n1) 315 k2=mod(n2+mod(i2-j2,n2),n2) 316 k3=mod(n3+mod(i3-j3,n3),n3) 317 ifft=k1+1+n1*(k2+n2*k3) 318 ! The signs of imaginary contributions have been checked 319 khxc(1,ipw1,1,ipw2,1)=kxcg(1,ifft) 320 khxc(2,ipw1,1,ipw2,1)=kxcg(2,ifft) 321 khxc(1,ipw2,1,ipw1,1)=kxcg(1,ifft) 322 khxc(2,ipw2,1,ipw1,1)=-kxcg(2,ifft) 323 end if 324 end do 325 end if 326 327 end if 328 end do 329 330 ABI_DEALLOCATE(wkxc) 331 ABI_DEALLOCATE(kxcg) 332 333 ! Endif option 2 334 end if 335 336 !Now, get the symmetric dielectric matrix 337 !Premultiplication by square root of minus susceptibility matrix 338 do ipw2=1,npwdiel 339 do ipw1=1,npwdiel 340 ar=0._dp ; ai=0._dp 341 do ii=1,npwdiel 342 ar=ar+sqrsus(1,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) & 343 & -sqrsus(2,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1) 344 ai=ai+sqrsus(2,ipw1,1,ii,1)*khxc(1,ii,1,ipw2,1) & 345 & +sqrsus(1,ipw1,1,ii,1)*khxc(2,ii,1,ipw2,1) 346 end do 347 work(1,ipw1,1,ipw2,1)=ar 348 work(2,ipw1,1,ipw2,1)=ai 349 end do 350 end do 351 !Postmultiplication by square root of minus susceptibility matrix 352 do ipw2=1,npwdiel 353 ! do ipw1=ipw2,npwdiel 354 do ipw1=1,npwdiel 355 ar=0._dp ; ai=0._dp 356 do ii=1,npwdiel 357 ar=ar+work(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 358 & -work(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 359 ai=ai+work(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 360 & +work(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 361 end do 362 symdielmat(1,ipw1,1,ipw2,1)=ar 363 symdielmat(2,ipw1,1,ipw2,1)=ai 364 ! if(ipw1/=ipw2)then 365 ! symdielmat(1,ipw2,1,ipw1,1)=ar 366 ! symdielmat(2,ipw2,1,ipw1,1)=-ai 367 ! end if 368 end do 369 ! Add the unity matrix 370 symdielmat(1,ipw2,1,ipw2,1)=1._dp+symdielmat(1,ipw2,1,ipw2,1) 371 end do 372 373 ABI_DEALLOCATE(khxc) 374 375 ABI_ALLOCATE(symh,(npwdiel*(npwdiel+1))) 376 ABI_ALLOCATE(symvec,(2,npwdiel,nspden,npwdiel,nspden)) 377 ABI_ALLOCATE(eig_sym,(npwdiel)) 378 379 !Store the symmetrized dielectric matrix in proper mode before calling zhpev 380 index=1 381 do ii=1,npwdiel 382 do jj=1,ii 383 symh(index )=symdielmat(1,jj,1,ii,1) 384 symh(index+1)=symdielmat(2,jj,1,ii,1) 385 index=index+2 386 end do 387 end do 388 389 ier=0 390 call ZHPEV ('V','U',npwdiel,symh,eig_sym,symvec,npwdiel,zhpev1,& 391 & zhpev2,ier) 392 393 if(prtvol>=10)then 394 write(message, '(a,a,a,5es12.4)' )ch10,& 395 & ' Five largest eigenvalues of the symmetrized dielectric matrix:',& 396 & ch10,eig_sym(npwdiel:npwdiel-4:-1) 397 call wrtout(ab_out,message,'COLL') 398 end if 399 400 write(message, '(a,a)' )ch10,& 401 & ' dieltcel : 15 largest eigenvalues of the symmetrized dielectric matrix' 402 call wrtout(std_out,message,'COLL') 403 write(message, '(a,5es12.5)' )' 1-5 :',eig_sym(npwdiel:npwdiel-4:-1) 404 call wrtout(std_out,message,'COLL') 405 write(message, '(a,5es12.5)' )' 6-10 :',eig_sym(npwdiel-5:npwdiel-9:-1) 406 call wrtout(std_out,message,'COLL') 407 write(message, '(a,5es12.5)' )' 11-15:',eig_sym(npwdiel-10:npwdiel-14:-1) 408 call wrtout(std_out,message,'COLL') 409 write(message, '(a,a)' )ch10,& 410 & ' dieltcel : 5 smallest eigenvalues of the symmetrized dielectric matrix' 411 call wrtout(std_out,message,'COLL') 412 write(message, '(a,5es12.5)' )' 1-5 :',eig_sym(1:5) 413 call wrtout(std_out,message,'COLL') 414 415 !Invert the hermitian dielectric matrix, 416 work(:,:,:,:,:)=0.0_dp 417 do ieig=1,npwdiel 418 eiginv=1.0_dp/eig_sym(ieig) 419 do ipw2=1,npwdiel 420 ! do ipw1=ipw2,npwdiel 421 do ipw1=1,npwdiel 422 work(1,ipw1,1,ipw2,1)=work(1,ipw1,1,ipw2,1)+& 423 & (symvec(1,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)+ & 424 & symvec(2,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv 425 work(2,ipw1,1,ipw2,1)=work(2,ipw1,1,ipw2,1)+& 426 & (symvec(2,ipw1,1,ieig,1)*symvec(1,ipw2,1,ieig,1)- & 427 & symvec(1,ipw1,1,ieig,1)*symvec(2,ipw2,1,ieig,1) ) * eiginv 428 end do 429 end do 430 end do 431 !if(npwdiel>1)then 432 !do ipw2=2,npwdiel 433 !do ipw1=1,ipw2-1 434 !work(1,ipw1,1,ipw2,1)= work(1,ipw2,1,ipw1,1) 435 !work(2,ipw1,1,ipw2,1)=-work(2,ipw2,1,ipw1,1) 436 !end do 437 !end do 438 !end if 439 440 ABI_DEALLOCATE(eig_sym) 441 ABI_DEALLOCATE(symh) 442 ABI_DEALLOCATE(symvec) 443 444 !DEBUG 445 !Checks whether the inverse of the symmetric dielectric matrix 446 !has been correctly generated 447 !do ipw1=1,npwdiel 448 !do ipw2=1,npwdiel 449 !elementr=0.0_dp 450 !elementi=0.0_dp 451 !do ipw3=1,npwdiel 452 !elementr=elementr+work(1,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1)& 453 !& -work(2,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1) 454 !elementi=elementi+work(1,ipw1,1,ipw3,1)*symdielmat(2,ipw3,1,ipw2,1)& 455 !& +work(2,ipw1,1,ipw3,1)*symdielmat(1,ipw3,1,ipw2,1) 456 !end do 457 !if(elementr**2+elementi**2 > 1.0d-12)then 458 !if( ipw1 /= ipw2 .or. & 459 !& ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then 460 !write(std_out,*)' dieltcel : the inversion procedure is not correct ' 461 !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2 462 !write(std_out,*)' elementr,elementi=',elementr,elementi 463 !stop 464 !end if 465 !end if 466 !end do 467 !end do 468 !write(std_out,*)'dieltcel : matrix has been inverted successfully ' 469 !ENDDEBUG 470 471 ABI_DEALLOCATE(symdielmat) 472 473 !Then get the inverse of the asymmetric 474 !dielectric matrix, as required for the preconditioning. 475 !Premultiplication by square root of minus susceptibility matrix 476 do ipw2=1,npwdiel 477 do ipw1=1,npwdiel 478 ar=0._dp ; ai=0._dp 479 do ii=1,npwdiel 480 ar=ar+invsqrsus(1,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) & 481 & -invsqrsus(2,ipw1,1,ii,1)*work(2,ii,1,ipw2,1) 482 ai=ai+invsqrsus(2,ipw1,1,ii,1)*work(1,ii,1,ipw2,1) & 483 & +invsqrsus(1,ipw1,1,ii,1)*work(2,ii,1,ipw2,1) 484 end do 485 work2(1,ipw1,1,ipw2,1)=ar 486 work2(2,ipw1,1,ipw2,1)=ai 487 end do 488 end do 489 !Postmultiplication by square root of minus susceptibility matrix 490 do ipw2=1,npwdiel 491 do ipw1=1,npwdiel 492 ar=0._dp ; ai=0._dp 493 do ii=1,npwdiel 494 ar=ar+work2(1,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 495 & -work2(2,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 496 ai=ai+work2(2,ipw1,1,ii,1)*sqrsus(1,ii,1,ipw2,1) & 497 & +work2(1,ipw1,1,ii,1)*sqrsus(2,ii,1,ipw2,1) 498 end do 499 dielinv(1,ipw1,1,ipw2,1)=ar 500 dielinv(2,ipw1,1,ipw2,1)=ai 501 end do 502 end do 503 504 ABI_DEALLOCATE(invsqrsus) 505 ABI_DEALLOCATE(sqrsus) 506 ABI_DEALLOCATE(work) 507 ABI_DEALLOCATE(work2) 508 ABI_DEALLOCATE(zhpev1) 509 ABI_DEALLOCATE(zhpev2) 510 511 call timab(96,2,tsec) 512 513 !DEBUG 514 !write(std_out,*)' dieltcel : exit, will stop' 515 !stop 516 !ENDDEBUG 517 518 end subroutine dieltcel