TABLE OF CONTENTS
ABINIT/suscep_stat [ Functions ]
NAME
suscep_stat
FUNCTION
Compute the susceptibility matrix from input wavefunctions, band occupations, and k point wts. Include the usual sum-over-state terms, but also the corrections due to the change of the Fermi level in the metallic case, as well as implicit sum over higher lying conduction states, thanks to the closure relation (referred to as an extrapolation).
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (XG,AR,MB) 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
atindx(natom)=index table for atoms (see scfcv.f) atindx1(natom)=index table for atoms, inverse of atindx (see scfcv.f) cg(2,mcg)=wf in G space cprj(natom,mcg*usecprj)= wave functions projected with non-local projectors: cprj_nk(i)=<p_i|Cnk> where p_i is a non-local projector. dielar(7)=input parameters for dielectric matrix and susceptibility: diecut,dielng,diemac,diemix,diegap,dielam,diemixmag. dimcprj(natom*usepaw)=array of dimensions of array cprj (ordered by atom-type) doccde(mband*nkpt*nsppol)=derivative of occupancies wrt the energy for each band and k point eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for the dielectric matrix gprimd(3,3)=dimensional reciprocal space primitive translations irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data istwfk(nkpt)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates. kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix. lmax_diel=1+max. value of l angular momentum used for dielectric matrix mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix mkmem=number of k points treated by this node mpi_enreg=information about MPI parallelization mpw=maximum allowed value for npw natom=number of atoms in cell nband(nkpt*nsppol)=number of bands to be included in summation at each k point for each spin channel neglect_pawhat=1 if PAW contribution from hat density (compensation charge) has to be neglected (to be used when only an estimation of suscep. matrix has to be evaluated, i.e. for SCF precondictioning) nfftdiel=number of fft grid points for the computation of the diel matrix ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix, see ~abinit/doc/input_variables/vargs.htm#ngfft nkpt=number of k points npwarr(nkpt)=number of planewaves and boundary planewaves at each k, for going from the WF sphere to the medium size FFT grid. npwdiel=third and fifth dimension of the susmat array. nspden=number of spin-density components nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetry elements in group (at least 1 for identity) ntypat=number of types of atoms in unit cell. occ(mband*nkpt*nsppol)= occupation numbers for each band (usually 2.0) at each k point occopt=option for occupancies pawang <type(pawang_type)>=paw angular mesh and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data phnonsdiel(2,nfftdiel**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases ph1ddiel(2,3*(2*mgfftdiel+1)*natom*usepaw)=one-dimensional structure factor information for the dielectric matrix rprimd(3,3)=dimensional real space primitive translations symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry matrices in real space (integers) tnons(3,nsym)=reduced nonsymmorphic translations (symrel and tnons are in terms of real space primitive translations) typat(natom)=type (integer) for each atom ucvol=unit cell volume (Bohr**3) unpaw=unit number for cprj PAW data (if used) usecprj= 1 if cprj array is stored in memory usepaw=flag for PAW usetimerev=1 if Time-Reversal symmetry has to be used when symmetrizing susceptibility wtk(nkpt)=k point weights (they sum to 1.0) ylmdiel(npwdiel,lmax_diel**2)= real spherical harmonics for each G and k point for the dielectric matrix
OUTPUT
susmat(2,npwdiel,nspden,npwdiel,nspden)= the susceptibility (or density-density response) matrix in reciprocal space
NOTES
Case of non-collinear magnetism: In principle, we should compute 16 susceptibility matrix: chi0-(s1,s2),(s3,s4) (where s1, s2, s3,and s4 are spin indexes)... But, for the time being, the susceptibility is only used to compute the dielectric matrix within RPA approximation; in this approximation, only four susceptibilities are non-zero: chi0-(s1,s1),(s3,s3). They are stored in susmat(:,ipw1,1:2,ipw2,1:2)
PARENTS
vtorho
CHILDREN
destroy_mpi_enreg,fftpac,init_distribfft_seq,initmpi_seq,pawcprj_alloc pawcprj_free,pawcprj_get,pawcprj_mpi_allgather,pawgylmg,ph1d3d sphereboundary,susk,suskmm,symg,symrhg,timab,xmpi_sum,zhpev
SOURCE
111 #if defined HAVE_CONFIG_H 112 #include "config.h" 113 #endif 114 115 #include "abi_common.h" 116 117 118 subroutine suscep_stat(atindx,atindx1,cg,cprj,dielar,dimcprj,doccde,& 119 & eigen,gbound_diel,gprimd,irrzondiel,istwfk,kg,& 120 & kg_diel,lmax_diel,& 121 & mband,mcg,mcprj,mgfftdiel,mkmem,mpi_enreg,mpw,natom,nband,& 122 & neglect_pawhat,nfftdiel,ngfftdiel,nkpt,npwarr,& 123 & npwdiel,nspden,nspinor,nsppol,nsym,ntypat,occ,occopt,& 124 & pawang,pawtab,phnonsdiel,ph1ddiel,rprimd,& 125 & susmat,symafm,symrel,tnons,typat,ucvol,unpaw,usecprj,usepaw,usetimerev,& 126 & wtk,ylmdiel) 127 128 use defs_basis 129 use defs_abitypes 130 use m_xmpi 131 use m_errors 132 use m_profiling_abi 133 134 use m_pawang, only : pawang_type 135 use m_pawtab, only : pawtab_type 136 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, & 137 & pawcprj_get, pawcprj_mpi_allgather, pawcprj_free 138 use m_mpinfo, only : destroy_mpi_enreg 139 140 !This section has been created automatically by the script Abilint (TD). 141 !Do not modify the following lines by hand. 142 #undef ABI_FUNC 143 #define ABI_FUNC 'suscep_stat' 144 use interfaces_18_timing 145 use interfaces_32_util 146 use interfaces_51_manage_mpi 147 use interfaces_52_fft_mpi_noabirule 148 use interfaces_53_ffts 149 use interfaces_56_recipspace 150 use interfaces_65_paw 151 use interfaces_67_common 152 use interfaces_77_suscep, except_this_one => suscep_stat 153 !End of the abilint section 154 155 implicit none 156 157 !Arguments ------------------------------------ 158 !scalars 159 integer,intent(in) :: lmax_diel,mband,mcg,mcprj,mgfftdiel,mkmem,mpw,natom,neglect_pawhat 160 integer,intent(in) :: nfftdiel,nkpt,npwdiel,nspden,nspinor,nsppol,nsym,ntypat,occopt 161 integer,intent(in) :: unpaw,usecprj,usepaw,usetimerev 162 real(dp),intent(in) :: ucvol 163 type(MPI_type),intent(in) :: mpi_enreg 164 type(pawang_type),intent(in) :: pawang 165 !arrays 166 integer,intent(in) :: atindx(natom),atindx1(natom),dimcprj(natom*usepaw) 167 integer,intent(in) :: gbound_diel(2*mgfftdiel+8,2) 168 !no_abirules 169 !nfftdiel**(1-1/nsym) is 1 if nsym==1, and nfftdiel otherwise 170 integer,intent(in) :: irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),& 171 & istwfk(nkpt) 172 integer,intent(in) :: kg(3,mpw*mkmem),kg_diel(3,npwdiel),& 173 & nband(nkpt*nsppol),ngfftdiel(18) 174 integer,intent(in) :: npwarr(nkpt),symafm(nsym),symrel(3,3,nsym),typat(ntypat) 175 real(dp),intent(in) :: cg(2,mcg),dielar(7) 176 real(dp),intent(in) :: doccde(mband*nkpt*nsppol),eigen(mband*nkpt*nsppol) 177 real(dp),intent(in) :: gprimd(3,3),occ(mband*nkpt*nsppol) 178 !nfftdiel**(1-1/nsym) is 1 if nsym==1, and nfftdiel otherwise 179 real(dp),intent(in) :: phnonsdiel(2,nfftdiel**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),& 180 & tnons(3,nsym),wtk(nkpt) 181 real(dp),intent(in) :: ph1ddiel(2,(3*(2*mgfftdiel+1)*natom)*usepaw),rprimd(3,3) 182 real(dp),intent(in) :: ylmdiel(npwdiel,lmax_diel**2) 183 real(dp),intent(out) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 184 type(pawcprj_type) :: cprj(natom,mcprj*usecprj) 185 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 186 187 !Local variables------------------------------- 188 !scalars 189 integer :: bdtot_index,diag,extrap,i1,i2,i3,iband,ibg,icg,ier,ierr 190 integer :: ifft,ii,ikg,ikpt,indx,iorder_cprj,ipw1,ipw2,isp,isp1,isp2 191 integer :: ispinor,istwf_k,isym,j1,j2,j3,jj,jsp,k1,k2,k3 192 integer :: my_nspinor,nband_k,nband_loc,ndiel1,ndiel2,ndiel3,ndiel4,ndiel5,ndiel6 193 integer :: nkpg_diel,npw_k,npwsp,nspden_eff,nspden_tmp,nsym1,nsym2 194 integer :: paral_kgb_diel,spaceComm,t1,t2,testocc 195 real(dp) :: ai,ai2,ar,ar2,diegap,dielam,emax,invnsym 196 real(dp) :: invnsym1,invnsym2,phi1,phi12,phi2,phr1,phr12 197 real(dp) :: phr2,sumdocc,weight 198 logical :: antiferro 199 character(len=500) :: message 200 type(MPI_type) :: mpi_enreg_diel 201 202 !arrays 203 integer,allocatable :: gbound(:,:),kg_k(:,:),sym_g(:,:) 204 integer,allocatable :: tmrev_g(:) 205 real(dp) :: kpt_diel(3,1),tsec(2) 206 real(dp),allocatable :: drhode(:,:,:),drhode_wk(:,:,:) 207 real(dp),allocatable :: eig_diel(:),gylmg_diel(:,:,:),kpg_dum(:,:) 208 real(dp),allocatable :: occ_deavg(:),ph3d_diel(:,:,:),phdiel(:,:,:) 209 real(dp),allocatable :: phkxred_diel(:,:),rhoextrap(:,:,:,:),rhoextrg(:,:) 210 real(dp),allocatable :: rhoextrr(:,:),sush(:),sussum(:),susvec(:,:,:) 211 real(dp),allocatable :: suswk(:,:,:,:),zhpev1(:,:),zhpev2(:) 212 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_loc(:,:) 213 214 ! ************************************************************************* 215 216 call timab(740,1,tsec) 217 call timab(741,1,tsec) 218 219 ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!") 220 221 222 !----- Initialisations ----------------------------------------------------------- 223 !--------------------------------------------------------------------------------- 224 225 if (usecprj==0.and.usepaw==1) then 226 write (message,'(3a)')& 227 & ' cprj datastructure must be allocated !',ch10,& 228 & ' Action: change pawusecp input keyword.' 229 MSG_ERROR(message) 230 end if 231 232 if (mpi_enreg%paral_spinor==1) then 233 message = ' not yet allowed for parallelization over spinors !' 234 MSG_ERROR(message) 235 end if 236 237 !Init mpicomm 238 if(mpi_enreg%paral_kgb==1) then 239 spaceComm=mpi_enreg%comm_kpt 240 else 241 spaceComm=mpi_enreg%comm_cell 242 end if 243 244 !The dielectric stuff is performed in sequential mode. 245 !Set mpi_enreg_diel accordingly 246 call initmpi_seq(mpi_enreg_diel) 247 call init_distribfft_seq(MPI_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all') 248 249 !testocc to be taken away 250 testocc=1 251 252 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor) 253 254 iorder_cprj=0 ! order for the cprj reading... 255 256 !Initialize some scalar quantities 257 antiferro=(nsppol==1.and.nspden==2) 258 nspden_eff=min(max(nsppol,nspden),2) ! Size for the computed part of susmat 259 bdtot_index=0 ; icg=0 ; ibg=0 260 ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3) 261 262 !ndiel4,ndiel5,ndiel6 are FFT dimensions, modified to avoid cache trashing 263 ndiel4=ngfftdiel(4) ; ndiel5=ngfftdiel(5) ; ndiel6=ngfftdiel(6) 264 diegap=dielar(5) ; dielam=dielar(6) 265 extrap=0 266 267 !If dielam is too small, there is no extrapolation. 268 if(dielam>1.0d-6)extrap=1 269 270 !Some stuff for symmetries 271 nsym1=sum(symafm,mask=symafm==1) 272 nsym2=nsym-nsym1 273 invnsym =one/dble(nsym) 274 invnsym1=one/dble(nsym1) 275 invnsym2=one 276 !FIXME: make sure this is consistent with following code 277 !div by 0 for several v5 tests 278 if (nsym2 > 0) invnsym2=one/dble(nsym2) 279 280 !Allocations 281 ABI_ALLOCATE(occ_deavg,(mband)) 282 if(occopt>=3) then 283 ABI_ALLOCATE(drhode,(2,npwdiel,nspden_eff)) 284 else 285 ABI_ALLOCATE(drhode,(0,0,0)) 286 end if 287 if(extrap==1) then 288 ABI_ALLOCATE(rhoextrap,(ndiel4,ndiel5,ndiel6,nspinor)) 289 else 290 ABI_ALLOCATE(rhoextrap,(0,0,0,0)) 291 end if 292 293 !zero the susceptibility matrix and other needed quantities 294 susmat(:,:,:,:,:)=zero 295 if(occopt>=3)then 296 drhode(:,:,:)=zero 297 sumdocc=zero 298 end if 299 300 !PAW additional initializations 301 if (usepaw==1) then 302 ABI_ALLOCATE(gylmg_diel,(npwdiel,lmax_diel**2,ntypat)) 303 ABI_ALLOCATE(ph3d_diel,(2,npwdiel,natom)) 304 if (neglect_pawhat==0) then 305 ABI_ALLOCATE(phkxred_diel,(2,natom)) 306 ABI_ALLOCATE(kpg_dum,(0,0)) 307 kpt_diel(1:3,1)=zero;phkxred_diel(1,:)=one;phkxred_diel(2,:)=zero;nkpg_diel=0 308 ! write(std_out,*) ' lmax_diel ', lmax_diel 309 call pawgylmg(gprimd,gylmg_diel,kg_diel,kpg_dum,kpt_diel,lmax_diel,nkpg_diel,npwdiel,ntypat,pawtab,ylmdiel) 310 call ph1d3d(1,natom,kg_diel,natom,natom,npwdiel,ndiel1,ndiel2,ndiel3,phkxred_diel,ph1ddiel,ph3d_diel) 311 ABI_DEALLOCATE(phkxred_diel) 312 ABI_DEALLOCATE(kpg_dum) 313 else 314 gylmg_diel=zero;ph3d_diel=one 315 end if 316 else 317 ABI_ALLOCATE(gylmg_diel,(0,0,0)) 318 ABI_ALLOCATE(ph3d_diel,(0,0,0)) 319 end if 320 321 call timab(741,2,tsec) 322 323 324 325 !--BIG loop over spins ------------------------------------------------------------ 326 !--------------------------------------------------------------------------------- 327 328 do isp=1,nsppol 329 ikg=0 330 331 if(extrap==1)rhoextrap(:,:,:,:)=zero 332 333 ! --BIG loop over k-points -------------------------------------------------------- 334 ! --------------------------------------------------------------------------------- 335 336 do ikpt=1,nkpt 337 338 nband_k=nband(ikpt+(isp-1)*nkpt) 339 istwf_k=istwfk(ikpt) 340 npw_k=npwarr(ikpt) 341 342 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isp,mpi_enreg%me_kpt)) then 343 bdtot_index=bdtot_index+nband_k 344 cycle 345 end if 346 347 call timab(742,1,tsec) 348 349 350 ABI_ALLOCATE(gbound,(2*mgfftdiel+8,2)) 351 ABI_ALLOCATE(kg_k,(3,npw_k)) 352 353 if (usepaw==1) then 354 ABI_DATATYPE_ALLOCATE(cprj_k,(natom,my_nspinor*nband_k)) 355 if (neglect_pawhat==0) then 356 call pawcprj_alloc(cprj_k,0,dimcprj) 357 if (mpi_enreg%nproc_band==1) then 358 call pawcprj_get(atindx1,cprj_k,cprj,natom,1,ibg,ikpt,iorder_cprj,isp,& 359 & mband,mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,unpaw,& 360 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 361 else 362 nband_loc=nband_k/mpi_enreg%nproc_band 363 ABI_DATATYPE_ALLOCATE(cprj_loc,(natom,my_nspinor*nband_loc)) 364 call pawcprj_alloc(cprj_loc,0,dimcprj) 365 call pawcprj_get(atindx1,cprj_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isp,& 366 & mband/mpi_enreg%nproc_band,mkmem,natom,nband_loc,nband_loc,my_nspinor,nsppol,unpaw,& 367 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 368 call pawcprj_mpi_allgather(cprj_loc,cprj_k,natom,my_nspinor*nband_loc,dimcprj,0,& 369 & mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.true.) 370 call pawcprj_free(cprj_loc) 371 ABI_DATATYPE_DEALLOCATE(cprj_loc) 372 end if 373 else 374 !call pawcprj_nullify(cprj_k) 375 end if 376 else 377 ABI_DATATYPE_ALLOCATE(cprj_k,(0,0)) 378 end if 379 380 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 381 call sphereboundary(gbound,istwf_k,kg_k,mgfftdiel,npw_k) 382 383 if(extrap==1)then 384 ! Compute inverse of average dielectric gap for each band 385 ! and multiply by occupation factor 386 emax=maxval(eigen(1+bdtot_index:nband_k+bdtot_index)) 387 do iband=1,nband_k 388 occ_deavg(iband)= occ(iband+bdtot_index)*dielam & 389 & / ( emax-eigen(iband+bdtot_index) + diegap ) 390 end do 391 else 392 occ_deavg(:)=zero 393 end if 394 395 call timab(742,2,tsec) 396 call timab(743,1,tsec) 397 398 ! Compute the contribution of each k-point to susmat, rhoextrap, drhode and sumdocc. 399 if(mpi_enreg%paral_kgb==1)then !Only this version is in parallel 400 ! Use either the simpler implementation 401 ! Should provide a test !!! 402 call susk(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,& 403 & gbound_diel,gylmg_diel,icg,ikpt,& 404 & isp,istwf_k,kg_diel,kg_k,lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,& 405 & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,& 406 & npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,& 407 & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,& 408 & susmat,typat,ucvol,usepaw,wtk) 409 else 410 ! Or the more sophisticated one, needed to save memory. 411 call suskmm(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,& 412 & gbound_diel,gylmg_diel,icg,ikpt,& 413 & isp,istwf_k,kg_diel,kg_k,lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,& 414 & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,& 415 & npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,paral_kgb_diel,& 416 & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,& 417 & susmat,typat,ucvol,usepaw,wtk) 418 end if 419 420 call timab(743,2,tsec) 421 422 ABI_DEALLOCATE(gbound) 423 ABI_DEALLOCATE(kg_k) 424 425 bdtot_index=bdtot_index+nband_k 426 427 if (mkmem/=0) then 428 ibg=ibg+my_nspinor*nband_k 429 icg=icg+my_nspinor*npw_k*nband_k 430 ikg=ikg+npw_k 431 end if 432 if (usepaw==1) then 433 if (neglect_pawhat==0) then 434 call pawcprj_free(cprj_k) 435 end if 436 end if 437 ABI_DATATYPE_DEALLOCATE(cprj_k) 438 439 ! End loop on ikpt: -------------------------------------------------------- 440 end do 441 442 ! Here include the contribution from the extrapolation to susmat, 443 ! diagonal part 444 if(extrap==1)then 445 446 call timab(744,1,tsec) 447 448 ! Transfer extrapolating density on augmented fft grid to 449 ! normal fft grid in real space. 450 ! Warning1 : if collinear magnetism, must treat only one spin at a time 451 ! Warning2 : if non-collinear magnetism, treat both spins 452 ! Warning3 : this is subtle for antiferro magnetism 453 nspden_tmp=1;if (antiferro) nspden_tmp=2 454 ABI_ALLOCATE(rhoextrr,(nfftdiel,nspden_tmp)) 455 ABI_ALLOCATE(rhoextrg,(2,nfftdiel)) 456 if (nspden==1.and.nspinor==2) rhoextrap(:,:,:,1)=rhoextrap(:,:,:,1)+rhoextrap(:,:,:,2) 457 458 do ispinor=1,min(nspinor,nspden) 459 jsp=isp+ispinor-1 460 461 call fftpac(1,mpi_enreg_diel,1,ndiel1,ndiel2,ndiel3,ndiel4,ndiel5,ndiel6,& 462 & ngfftdiel,rhoextrr(:,1),rhoextrap(:,:,:,ispinor),1) 463 464 ! Generate the density in reciprocal space, and symmetrize it 465 ! (note symrhg also make the reverse FFT, to get symmetrized density; 466 ! this is useless here, and should be made an option) 467 call symrhg(1,gprimd,irrzondiel,mpi_enreg_diel,nfftdiel,nfftdiel,ngfftdiel,& 468 & nspden_tmp,1,nsym,paral_kgb_diel,phnonsdiel,rhoextrg,rhoextrr,rprimd,symafm,symrel) 469 470 do ipw2=1,npwdiel 471 j1=kg_diel(1,ipw2) ; j2=kg_diel(2,ipw2) ; j3=kg_diel(3,ipw2) 472 ! static: Only fills lower half of the matrix (here, the susceptibility matrix) 473 ! dynamical: fill all, will not affect susopt==2 for which extrap==0 474 do ipw1=1,npwdiel 475 i1=kg_diel(1,ipw1) ; i2=kg_diel(2,ipw1) ; i3=kg_diel(3,ipw1) 476 ! NOTE that there is a FFT folding (superposition) bias here 477 ! Should use kgindex, in the same spirit as in prcref 478 k1=i1-j1; k1=modulo(k1,ndiel1) 479 k2=i2-j2; k2=modulo(k2,ndiel2) 480 k3=i3-j3; k3=modulo(k3,ndiel3) 481 ifft=k1+1+ndiel1*(k2+ndiel2*k3) 482 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+rhoextrg(1,ifft) 483 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+rhoextrg(2,ifft) 484 end do 485 end do 486 487 end do 488 ABI_DEALLOCATE(rhoextrg) 489 ABI_DEALLOCATE(rhoextrr) 490 491 call timab(744,2,tsec) 492 493 end if 494 495 ! End loop over spins --------------------------------------------------------- 496 end do 497 498 ABI_DEALLOCATE(occ_deavg) 499 ABI_DEALLOCATE(rhoextrap) 500 ABI_DEALLOCATE(gylmg_diel) 501 ABI_DEALLOCATE(ph3d_diel) 502 !end if 503 504 call destroy_mpi_enreg(mpi_enreg_diel) 505 506 !-- Stuff for parallelism -------------------------------------------------------- 507 !--------------------------------------------------------------------------------- 508 509 if(xmpi_paral==1)then 510 call timab(746,1,tsec) 511 ABI_ALLOCATE(sussum,(2*npwdiel*nspden*npwdiel*nspden)) 512 ! Recreate full susmat on all proc. 513 ! This should be coded more efficiently, 514 ! since half of the matrix is still empty, and 515 ! it is spin-diagonal. 516 sussum(:)=reshape(susmat(:,:,:,:,:),(/2*npwdiel*nspden*npwdiel*nspden/)) 517 call xmpi_sum(sussum,spaceComm,ierr) 518 susmat(:,:,:,:,:)=reshape(sussum(:),(/2,npwdiel,nspden,npwdiel,nspden/)) 519 ABI_DEALLOCATE(sussum) 520 ! Recreate full drhode on all proc. 521 if(occopt>=3 .and. testocc==1)then 522 call xmpi_sum(drhode,spaceComm,ierr) 523 ! Should use only one mpi-allreduce call instead of the three 524 call xmpi_sum(sumdocc,spaceComm,ierr) 525 end if 526 call timab(746,2,tsec) 527 end if 528 529 !-- Apply spatial hermitian/symmetries on spin-diagonal susceptibility matrix ---- 530 !--------------------------------------------------------------------------------- 531 532 call timab(747,1,tsec) 533 534 !If antiferro magnetism, has to divide (spin-diagonal) susceptibility by 2 (due to dble occupations) 535 if (antiferro) then 536 do ipw2=1,npwdiel 537 do ipw1=ipw2,npwdiel 538 susmat(:,ipw1,1,ipw2,1)=half*susmat(:,ipw1,1,ipw2,1) 539 end do 540 end do 541 end if 542 543 !Generate upper half of the spin-diagonal matrix (still the susceptibility matrix) 544 do isp=1,nspden_eff 545 do ipw2=2,npwdiel 546 do ipw1=1,ipw2-1 547 susmat(1,ipw1,isp,ipw2,isp)= susmat(1,ipw2,isp,ipw1,isp) 548 susmat(2,ipw1,isp,ipw2,isp)=-susmat(2,ipw2,isp,ipw1,isp) 549 end do 550 end do 551 end do 552 553 !Compute symmetric of G-vectors and eventual phases 554 !(either time-reversal or spatial symmetries) 555 ABI_ALLOCATE(tmrev_g,(npwdiel)) 556 ABI_ALLOCATE(sym_g,(npwdiel,nsym)) 557 ABI_ALLOCATE(phdiel,(2,npwdiel,nsym)) 558 call symg(kg_diel,npwdiel,nsym,phdiel,sym_g,symrel,tmrev_g,tnons) 559 560 !Impose spatial symmetries to the spin-diagonal susceptibility matrix 561 ABI_ALLOCATE(suswk,(2,npwdiel,npwdiel,nspden_eff)) 562 do isp=1,nspden_eff 563 suswk(:,:,:,isp)=susmat(:,:,isp,:,isp) ! Temporary storage 564 end do 565 566 do isp=1,nspden_eff 567 jsp=min(3-isp,nsppol) 568 do ipw2=1,npwdiel 569 do ipw1=1,npwdiel 570 ar=suswk(1,ipw1,ipw2,isp) 571 ai=suswk(2,ipw1,ipw2,isp) 572 ar2=zero;ai2=zero 573 if(nsym>1)then 574 do isym=2,nsym 575 t1=sym_g(ipw1,isym) ; t2=sym_g(ipw2,isym) 576 ! Not all symmetries are non-symmorphic. Should save time here ... 577 phr1=phdiel(1,ipw1,isym) ; phi1=phdiel(2,ipw1,isym) 578 phr2=phdiel(1,ipw2,isym) ; phi2=phdiel(2,ipw2,isym) 579 phr12= phr1*phr2+phi1*phi2 ; phi12=phi1*phr2-phr1*phi2 580 if (symafm(isym)==1) then 581 ar=ar+suswk(1,t1,t2,isp)*phr12-suswk(2,t1,t2,isp)*phi12 582 ai=ai+suswk(2,t1,t2,isp)*phr12+suswk(1,t1,t2,isp)*phi12 583 else 584 ar2=ar2+suswk(1,t1,t2,jsp)*phr12-suswk(2,t1,t2,jsp)*phi12 585 ai2=ai2+suswk(2,t1,t2,jsp)*phr12+suswk(1,t1,t2,jsp)*phi12 586 end if 587 end do 588 end if 589 if (antiferro) then 590 susmat(1,ipw1,1,ipw2,1)=ar*invnsym1 591 susmat(2,ipw1,1,ipw2,1)=ai*invnsym1 592 susmat(1,ipw1,2,ipw2,2)=ar2*invnsym2 593 susmat(2,ipw1,2,ipw2,2)=ai2*invnsym2 594 else 595 susmat(1,ipw1,isp,ipw2,isp)=(ar+ar2)*invnsym 596 susmat(2,ipw1,isp,ipw2,isp)=(ai+ai2)*invnsym 597 end if 598 end do 599 end do 600 end do 601 ABI_DEALLOCATE(suswk) 602 603 604 !-- Add contribibution to susceptibility due to change of Fermi level ----------- 605 !--------------------------------------------------------------------------------- 606 607 if (occopt>=3.and.testocc==1) then 608 609 ! Impose spatial symmetries to drhode 610 ABI_ALLOCATE(drhode_wk,(2,npwdiel,nspden_eff)) 611 do isp=1,nspden_eff 612 jsp=min(3-isp,nsppol) 613 do ipw1=1,npwdiel 614 ar=drhode(1,ipw1,isp) 615 ai=drhode(2,ipw1,isp) 616 ar2=zero;ai2=zero 617 if (nsym>1) then 618 do isym=2,nsym 619 t1=sym_g(ipw1,isym) 620 ! Not all symmetries are non-symmorphic. Should save time here ... 621 phr1=phdiel(1,ipw1,isym);phi1=phdiel(2,ipw1,isym) 622 if (symafm(isym)==1) then 623 ar=ar+drhode(1,t1,isp)*phr1-drhode(2,t1,isp)*phi1 624 ai=ai+drhode(2,t1,isp)*phr1+drhode(1,t1,isp)*phi1 625 else 626 ar2=ar2+drhode(1,t1,jsp)*phr1-drhode(2,t1,jsp)*phi1 627 ai2=ai2+drhode(2,t1,jsp)*phr1+drhode(1,t1,jsp)*phi1 628 end if 629 end do 630 end if 631 if (antiferro) then ! 1/2 factor due to (dble) occupations 632 drhode_wk(1,ipw1,1)=half*ar*invnsym1 633 drhode_wk(2,ipw1,1)=half*ai*invnsym1 634 drhode_wk(1,ipw1,2)=half*ar2*invnsym2 635 drhode_wk(2,ipw1,2)=half*ai2*invnsym2 636 else 637 drhode_wk(1,ipw1,isp)=(ar+ar2)*invnsym 638 drhode_wk(2,ipw1,isp)=(ai+ai2)*invnsym 639 end if 640 end do 641 end do 642 643 ! Add contribution to non-diagonal susceptibility 644 ! Presently fills complete susceptibility matrix, not only lower half 645 weight=one/sumdocc 646 do isp2=1,nspden_eff 647 do ipw2=1,npwdiel 648 do isp1=1,nspden_eff 649 do ipw1=1,npwdiel 650 susmat(1,ipw1,isp1,ipw2,isp2)=susmat(1,ipw1,isp1,ipw2,isp2)- & 651 & weight*( drhode_wk(1,ipw1,isp1)*drhode_wk(1,ipw2,isp2) & 652 & +drhode_wk(2,ipw1,isp1)*drhode_wk(2,ipw2,isp2) ) 653 susmat(2,ipw1,isp1,ipw2,isp2)=susmat(2,ipw1,isp1,ipw2,isp2)- & 654 & weight*( drhode_wk(2,ipw1,isp1)*drhode_wk(1,ipw2,isp2) & 655 & -drhode_wk(1,ipw1,isp1)*drhode_wk(2,ipw2,isp2) ) 656 end do 657 end do 658 end do 659 end do 660 ABI_DEALLOCATE(drhode_wk) 661 662 end if 663 !if (occopt>=3) then 664 ABI_DEALLOCATE(drhode) 665 !end if 666 667 668 !--- Impose the time-reversal symmetry to the susceptibility matrix -------------- 669 !--------------------------------------------------------------------------------- 670 671 if (usetimerev==1) then 672 ABI_ALLOCATE(suswk,(2,npwdiel,npwdiel,1)) 673 674 ! Impose the time-reversal symmetry to the spin-diagonal susceptibility matrix 675 do isp=1,nspden_eff 676 suswk(:,:,:,1)=susmat(:,:,isp,:,isp) ! Temporary storage 677 do ipw2=1,npwdiel 678 t2=tmrev_g(ipw2) 679 do ipw1=1,npwdiel 680 t1=tmrev_g(ipw1) 681 susmat(1,ipw1,isp,ipw2,isp)=half*(suswk(1,ipw1,ipw2,1)+suswk(1,t1,t2,1)) 682 susmat(2,ipw1,isp,ipw2,isp)=half*(suswk(2,ipw1,ipw2,1)-suswk(2,t1,t2,1)) 683 end do 684 end do 685 end do 686 687 ! Impose the time-reversal symmetry to the off-diagonal susceptibility matrix 688 if (nspden_eff/=1.and.occopt>=3.and.testocc==1) then 689 suswk(:,:,:,1)=susmat(:,:,1,:,2) ! Temporary storage 690 do ipw2=1,npwdiel 691 t2=tmrev_g(ipw2) 692 do ipw1=1,npwdiel 693 t1=tmrev_g(ipw1) 694 ar=half*(suswk(1,ipw1,ipw2,1)+suswk(1,t1,t2,1)) 695 ai=half*(suswk(2,ipw1,ipw2,1)-suswk(2,t1,t2,1)) 696 susmat(1,ipw1,1,ipw2,2)= ar 697 susmat(2,ipw1,1,ipw2,2)= ai 698 susmat(1,ipw1,2,ipw2,1)= ar 699 susmat(2,ipw1,2,ipw2,1)=-ai 700 end do 701 end do 702 end if 703 ABI_DEALLOCATE(suswk) 704 end if 705 706 ABI_DEALLOCATE(phdiel) 707 ABI_DEALLOCATE(sym_g) 708 ABI_DEALLOCATE(tmrev_g) 709 710 711 !-- The full susceptibility matrix is computed ----------------------------------- 712 !-- Now, eventually diagonalize it and stop -------------------------------------- 713 !--------------------------------------------------------------------------------- 714 715 !Must turn on this flag to make the diagonalisation 716 diag=0 717 if(diag==1)then 718 719 npwsp=npwdiel*nspden_eff 720 ABI_ALLOCATE(sush,(npwsp*(npwsp+1))) 721 ABI_ALLOCATE(susvec,(2,npwsp,npwsp)) 722 ABI_ALLOCATE(eig_diel,(npwsp)) 723 ABI_ALLOCATE(zhpev1,(2,2*npwsp-1)) 724 ABI_ALLOCATE(zhpev2,(3*npwsp-2)) 725 ier=0 726 727 ! Store the susceptibility matrix in proper mode before calling zhpev 728 indx=1 729 do ii=1,npwdiel 730 do jj=1,ii 731 sush(indx )=susmat(1,jj,1,ii,1) 732 sush(indx+1)=susmat(2,jj,1,ii,1) 733 indx=indx+2 734 end do 735 end do 736 737 ! If spin-polarized, need to store other parts of the matrix 738 if(nspden_eff/=1)then 739 do ii=1,npwdiel 740 ! Here, spin-flip contribution 741 do jj=1,npwdiel 742 sush(indx )=susmat(1,jj,1,ii,2) 743 sush(indx+1)=susmat(2,jj,1,ii,2) 744 indx=indx+2 745 end do 746 ! Here spin down-spin down upper matrix 747 do jj=1,ii 748 sush(indx )=susmat(1,jj,2,ii,2) 749 sush(indx+1)=susmat(2,jj,2,ii,2) 750 indx=indx+2 751 end do 752 end do 753 end if 754 755 call ZHPEV ('V','U',npwsp,sush,eig_diel,susvec,npwsp,zhpev1,& 756 & zhpev2,ier) 757 758 write(std_out,*)' suscep_stat : print eigenvalues of the susceptibility matrix' 759 do ii=1,npwsp 760 write(std_out,'(i5,es16.6)' )ii,eig_diel(ii) 761 end do 762 763 ABI_DEALLOCATE(sush) 764 ABI_DEALLOCATE(susvec) 765 ABI_DEALLOCATE(eig_diel) 766 ABI_DEALLOCATE(zhpev1) 767 ABI_DEALLOCATE(zhpev2) 768 MSG_ERROR("Stopping here!") 769 end if 770 771 call timab(747,2,tsec) 772 call timab(740,2,tsec) 773 774 end subroutine suscep_stat