TABLE OF CONTENTS
ABINIT/dielmt [ Functions ]
NAME
dielmt
FUNCTION
Compute dielectric matrix 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. npwdiel=size of the dielinv and susmat arrays. nspden=number of spin-density components occopt=option for occupancies 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+4)/3,npwdiel,(nspden+4)/3)=inverse of the (non-hermitian)
TC dielectric matrix in reciprocal space.
NOTES
Warning : will not work in the spin-polarized, metallic case. Output (not cleaned) !!! Spin behaviour is not obvious !!!
TODO
Write equation below (hermitian matrix)
PARENTS
prcref,prcref_PMA
CHILDREN
timab,wrtout,zhpev
SOURCE
49 #if defined HAVE_CONFIG_H 50 #include "config.h" 51 #endif 52 53 #include "abi_common.h" 54 55 subroutine dielmt(dielinv,gmet,kg_diel,& 56 & npwdiel,nspden,occopt,prtvol,susmat) 57 58 use m_profiling_abi 59 60 use defs_basis 61 62 !This section has been created automatically by the script Abilint (TD). 63 !Do not modify the following lines by hand. 64 #undef ABI_FUNC 65 #define ABI_FUNC 'dielmt' 66 use interfaces_14_hidewrite 67 use interfaces_18_timing 68 !End of the abilint section 69 70 implicit none 71 72 !Arguments ------------------------------------ 73 !scalars 74 integer,intent(in) :: npwdiel,nspden,occopt,prtvol 75 !arrays 76 integer,intent(in) :: kg_diel(3,npwdiel) 77 real(dp),intent(in) :: gmet(3,3) 78 real(dp),intent(in) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 79 real(dp),intent(out) :: dielinv(2,npwdiel,nspden,npwdiel,nspden) 80 81 !Local variables------------------------------- 82 !scalars 83 integer :: ieig,ier,ii,index,ipw,ipw1,ipw2,isp,jj,npwsp 84 real(dp) :: ai1,ai2,ar1,ar2,eiginv,gfact,gfactinv,gred1,gred2,gred3,gsquar 85 real(dp) :: tpisq 86 character(len=500) :: message 87 !arrays 88 real(dp) :: tsec(2) 89 real(dp),allocatable :: dielh(:),dielmat(:,:,:,:,:),dielvec(:,:,:) 90 real(dp),allocatable :: eig_diel(:),zhpev1(:,:),zhpev2(:) 91 !no_abirules 92 !integer :: ipw3 93 !real(dp) :: elementi,elementr 94 95 ! ************************************************************************* 96 97 !DEBUG 98 !write(std_out,*)' dielmt : enter ' 99 !if(.true.)stop 100 !ENDDEBUG 101 102 !tpisq is (2 Pi) **2: 103 tpisq=(two_pi)**2 104 105 call timab(90,1,tsec) 106 107 !-Compute now the hermitian dielectric matrix------------------------------ 108 !Following remarks are only valid within RPA approximation (Kxc=0): 109 110 !for the spin-unpolarized case, 1 - 4pi (1/G) chi0(G,Gp) (1/Gp) 111 112 !for the spin-polarized case, 113 !( 1 0 ) - 4pi ( 1/G 1/G ) ( chi0 upup chi0 updn ) ( 1/Gp 1/Gp ) 114 !( 0 1 ) ( 1/G 1/G ) ( chi0 dnup chi0 dndn ) ( 1/Gp 1/Gp ) 115 !which is equal to 116 !( 1 0 ) - 4pi (1/G 0 ) (chi0 upup+dndn+updn+dnup chi0 upup+dndn+updn+dnup) (1/Gp 0 ) 117 !( 0 1 ) ( 0 1/G) (chi0 upup+dndn+updn+dnup chi0 upup+dndn+updn+dnup) ( 0 1/Gp) 118 !So, if spin-polarized, sum all spin contributions 119 !Note: chi0 updn = chi0 dnup = zero for non-metallic systems 120 121 !In the case of non-collinear magnetism, within RPA, this is the same because: 122 !chi0_(s1,s2),(s3,s4) = delta_s1,s2 * delta_s3,s4 * chi0_(s1,s1),(s3,s3) 123 !Only chi_upup,upup, chi_dndn,dndn, chi_upup,dndn and chi_dndn,upup 124 !have to be taken into account (stored, susmat(:,ipw1,1:2,ipw2,1:2) 125 126 ABI_ALLOCATE(dielmat,(2,npwdiel,min(nspden,2),npwdiel,min(nspden,2))) 127 128 if(nspden/=1)then 129 if (occopt<3) then 130 do ipw2=1,npwdiel 131 do ipw1=1,npwdiel 132 dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)+susmat(1,ipw1,2,ipw2,2) 133 dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)+susmat(2,ipw1,2,ipw2,2) 134 end do 135 end do 136 else 137 do ipw2=1,npwdiel 138 do ipw1=1,npwdiel 139 dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1)+susmat(1,ipw1,2,ipw2,2)+susmat(1,ipw1,1,ipw2,2)+susmat(1,ipw1,2,ipw2,1) 140 dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1)+susmat(2,ipw1,2,ipw2,2)+susmat(2,ipw1,1,ipw2,2)+susmat(2,ipw1,2,ipw2,1) 141 end do 142 end do 143 end if 144 else 145 do ipw2=1,npwdiel 146 do ipw1=1,npwdiel 147 dielmat(1,ipw1,1,ipw2,1)=susmat(1,ipw1,1,ipw2,1) 148 dielmat(2,ipw1,1,ipw2,1)=susmat(2,ipw1,1,ipw2,1) 149 end do 150 end do 151 end if 152 !Compute 1/G factors and include them in the dielectric matrix 153 do ipw1=1,npwdiel 154 gred1=dble(kg_diel(1,ipw1)) 155 gred2=dble(kg_diel(2,ipw1)) 156 gred3=dble(kg_diel(3,ipw1)) 157 gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 & 158 & +two*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 + & 159 & gmet(2,3)*gred2*gred3) ) 160 ! Distinguish G=0 from other elements 161 if(gsquar>tol12)then 162 ! !$ gfact=\sqrt (4.0_dp \pi/gsquar/dble(nspden))$ 163 gfact=sqrt(four_pi/gsquar) 164 do ipw2=1,npwdiel 165 ! Must multiply both rows and columns, and also changes the sign 166 dielmat(1,ipw2,1,ipw1,1)=-dielmat(1,ipw2,1,ipw1,1)*gfact 167 dielmat(2,ipw2,1,ipw1,1)=-dielmat(2,ipw2,1,ipw1,1)*gfact 168 dielmat(1,ipw1,1,ipw2,1)= dielmat(1,ipw1,1,ipw2,1)*gfact 169 dielmat(2,ipw1,1,ipw2,1)= dielmat(2,ipw1,1,ipw2,1)*gfact 170 end do 171 else 172 ! Zero the G=0 elements, head and wings 173 do ipw2=1,npwdiel 174 dielmat(1,ipw2,1,ipw1,1)=zero 175 dielmat(2,ipw2,1,ipw1,1)=zero 176 dielmat(1,ipw1,1,ipw2,1)=zero 177 dielmat(2,ipw1,1,ipw2,1)=zero 178 end do 179 end if 180 end do 181 182 !Complete the matrix in the spin-polarized case 183 !should this be nspden==2?? 184 if(nspden/=1)then 185 do ipw1=1,npwdiel 186 do ipw2=1,npwdiel 187 dielmat(1,ipw1,1,ipw2,2)=dielmat(1,ipw1,1,ipw2,1) 188 dielmat(2,ipw1,1,ipw2,2)=dielmat(2,ipw1,1,ipw2,1) 189 dielmat(1,ipw1,2,ipw2,1)=dielmat(1,ipw1,1,ipw2,1) 190 dielmat(2,ipw1,2,ipw2,1)=dielmat(2,ipw1,1,ipw2,1) 191 dielmat(1,ipw1,2,ipw2,2)=dielmat(1,ipw1,1,ipw2,1) 192 dielmat(2,ipw1,2,ipw2,2)=dielmat(2,ipw1,1,ipw2,1) 193 end do 194 end do 195 end if 196 197 !DEBUG 198 !write(std_out,*)' dielmt : make dielmat equal to identity matrix ' 199 !do ipw1=1,npwdiel 200 !do ipw2=1,npwdiel 201 !dielmat(1,ipw1,1,ipw2,1)=0.0_dp 202 !dielmat(2,ipw1,1,ipw2,1)=0.0_dp 203 !end do 204 !end do 205 !ENDDEBUG 206 207 !Add the diagonal part 208 do isp=1,min(nspden,2) 209 do ipw=1,npwdiel 210 dielmat(1,ipw,isp,ipw,isp)=one+dielmat(1,ipw,isp,ipw,isp) 211 end do 212 end do 213 214 !-The hermitian dielectric matrix is computed ------------------------------ 215 !-Now, diagonalize it ------------------------------------------------------ 216 217 !In RPA, everything is projected on the spin-symmetrized 218 !space. This was coded here (for the time being). 219 220 !Diagonalize the hermitian dielectric matrix 221 222 !npwsp=npwdiel*nspden 223 npwsp=npwdiel 224 225 ABI_ALLOCATE(dielh,(npwsp*(npwsp+1))) 226 ABI_ALLOCATE(dielvec,(2,npwsp,npwsp)) 227 ABI_ALLOCATE(eig_diel,(npwsp)) 228 ABI_ALLOCATE(zhpev1,(2,2*npwsp-1)) 229 ABI_ALLOCATE(zhpev2,(3*npwsp-2)) 230 ier=0 231 !Store the dielectric matrix in proper mode before calling zhpev 232 index=1 233 do ii=1,npwdiel 234 do jj=1,ii 235 dielh(index )=dielmat(1,jj,1,ii,1) 236 dielh(index+1)=dielmat(2,jj,1,ii,1) 237 index=index+2 238 end do 239 end do 240 !If spin-polarized and non RPA, need to store other parts of the matrix 241 !if(nspden/=1)then 242 !do ii=1,npwdiel 243 !Here, spin-flip contribution 244 !do jj=1,npwdiel 245 !dielh(index )=dielmat(1,jj,1,ii,2) 246 !dielh(index+1)=dielmat(2,jj,1,ii,2) 247 !index=index+2 248 !end do 249 !Here spin down-spin down upper matrix 250 !do jj=1,ii 251 !dielh(index )=dielmat(1,jj,2,ii,2) 252 !dielh(index+1)=dielmat(2,jj,2,ii,2) 253 !index=index+2 254 !end do 255 !end do 256 !end if 257 258 call ZHPEV ('V','U',npwsp,dielh,eig_diel,dielvec,npwdiel,zhpev1,& 259 & zhpev2,ier) 260 ABI_DEALLOCATE(zhpev1) 261 ABI_DEALLOCATE(zhpev2) 262 263 if(prtvol>=10)then 264 write(message, '(a,a,a,5es12.4)' )ch10,& 265 & ' Five largest eigenvalues of the hermitian RPA dielectric matrix:',& 266 & ch10,eig_diel(npwdiel:npwdiel-4:-1) 267 call wrtout(ab_out,message,'COLL') 268 end if 269 270 write(message, '(a,a)' )ch10,& 271 & ' dielmt : 15 largest eigenvalues of the hermitian RPA dielectric matrix' 272 call wrtout(std_out,message,'COLL') 273 write(message, '(a,5es12.5)' )' 1-5 :',eig_diel(npwdiel:npwdiel-4:-1) 274 call wrtout(std_out,message,'COLL') 275 write(message, '(a,5es12.5)' )' 6-10 :',eig_diel(npwdiel-5:npwdiel-9:-1) 276 call wrtout(std_out,message,'COLL') 277 write(message, '(a,5es12.5)' )' 11-15:',eig_diel(npwdiel-10:npwdiel-14:-1) 278 call wrtout(std_out,message,'COLL') 279 write(message, '(a,a)' )ch10,& 280 & ' dielmt : 5 smallest eigenvalues of the hermitian RPA dielectric matrix' 281 call wrtout(std_out,message,'COLL') 282 write(message, '(a,5es12.5)' )' 1-5 :',eig_diel(1:5) 283 call wrtout(std_out,message,'COLL') 284 285 !Invert the hermitian dielectric matrix, 286 !Should use a BLAS call ! 287 do ipw2=1,npwdiel 288 do ipw1=ipw2,npwdiel 289 dielinv(1,ipw1,1,ipw2,1)=zero 290 dielinv(2,ipw1,1,ipw2,1)=zero 291 end do 292 end do 293 do ieig=1,npwdiel 294 eiginv=one/eig_diel(ieig) 295 do ipw2=1,npwdiel 296 do ipw1=ipw2,npwdiel 297 ar1=dielvec(1,ipw1,ieig) 298 ai1=dielvec(2,ipw1,ieig) 299 ar2=dielvec(1,ipw2,ieig) 300 ai2=dielvec(2,ipw2,ieig) 301 dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)+& 302 & (ar1*ar2+ai1*ai2)*eiginv 303 dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)+& 304 & (ai1*ar2-ar1*ai2)*eiginv 305 end do 306 end do 307 end do 308 do ipw2=1,npwdiel-1 309 do ipw1=ipw2+1,npwdiel 310 dielinv(1,ipw2,1,ipw1,1)= dielinv(1,ipw1,1,ipw2,1) 311 dielinv(2,ipw2,1,ipw1,1)=-dielinv(2,ipw1,1,ipw2,1) 312 end do 313 end do 314 315 ABI_DEALLOCATE(dielh) 316 ABI_DEALLOCATE(dielvec) 317 ABI_DEALLOCATE(eig_diel) 318 319 !DEBUG 320 !Checks whether the inverse of the hermitian dielectric matrix 321 !has been correctly generated 322 !do ipw1=1,npwdiel 323 !do ipw2=1,npwdiel 324 !elementr=0.0_dp 325 !elementi=0.0_dp 326 !do ipw3=1,npwdiel 327 !elementr=elementr+dielinv(1,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1)& 328 !& -dielinv(2,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1) 329 !elementi=elementi+dielinv(1,ipw1,1,ipw3,1)*dielmat(2,ipw3,1,ipw2,1)& 330 !& +dielinv(2,ipw1,1,ipw3,1)*dielmat(1,ipw3,1,ipw2,1) 331 !end do 332 !if(elementr**2+elementi**2 > 1.0d-12)then 333 !if( ipw1 /= ipw2 .or. & 334 !& ( abs(elementr-1.0_dp)>1.0d-6 .or. abs(elementi)>1.0d-6 ))then 335 !write(std_out,*)' dielmt : the inversion procedure is not correct ' 336 !write(std_out,*)' ipw1, ipw2 =',ipw1,ipw2 337 !write(std_out,*)' elementr,elementi=',elementr,elementi 338 !stop 339 !end if 340 !end if 341 !end do 342 !end do 343 !write(std_out,*)'dielmt : matrix has been inverted successfully ' 344 !stop 345 !ENDDEBUG 346 347 !Then get the inverse of the asymmetric 348 !dielectric matrix, as required for the preconditioning. 349 350 !Inverse of the dielectric matrix : ( 1 - 4pi (1/G^2) chi0(G,Gp) )^(-1) 351 !In dielinv there is now (1 - 4pi (1/G) chi0(G,Gp) (1/Gp) )^(-1) 352 !So, evaluate dielinv_after(G,Gp) = 353 !(4pi/G^2)^(1/2) dielinv_before(G,Gp) (4pi/Gp^2)^(-1/2) 354 !In RPA, can focus on the spin-averaged quantities 355 do ipw1=1,npwdiel 356 gred1=dble(kg_diel(1,ipw1)) 357 gred2=dble(kg_diel(2,ipw1)) 358 gred3=dble(kg_diel(3,ipw1)) 359 gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 & 360 & +two*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 + & 361 & gmet(2,3)*gred2*gred3) ) 362 ! Distinguish G=0 from other elements 363 if(gsquar>tol12)then 364 gfact=sqrt(four_pi/gsquar) 365 gfactinv=one/gfact 366 do ipw2=1,npwdiel 367 ! Must multiply both rows and columns 368 dielinv(1,ipw2,1,ipw1,1)=dielinv(1,ipw2,1,ipw1,1)*gfactinv 369 dielinv(2,ipw2,1,ipw1,1)=dielinv(2,ipw2,1,ipw1,1)*gfactinv 370 dielinv(1,ipw1,1,ipw2,1)=dielinv(1,ipw1,1,ipw2,1)*gfact 371 dielinv(2,ipw1,1,ipw2,1)=dielinv(2,ipw1,1,ipw2,1)*gfact 372 end do 373 else 374 ! Zero the G=0 elements, head 375 do ipw2=1,npwdiel 376 if (ipw2/=ipw1) dielinv(1:2,ipw1,1,ipw2,1)=zero 377 end do 378 end if 379 end do 380 381 ABI_DEALLOCATE(dielmat) 382 383 call timab(90,2,tsec) 384 385 end subroutine dielmt