Line data Source code
1 : !!****m* ABINIT/m_suscep_stat
2 : !! NAME
3 : !! m_suscep_stat
4 : !!
5 : !! FUNCTION
6 : !! Compute the susceptibility matrix
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (XG, AR, MB)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : MODULE m_suscep_stat
23 :
24 : use defs_basis
25 : use m_xmpi
26 : use m_errors
27 : use m_abicore
28 : use m_distribfft
29 :
30 : use defs_abitypes, only : MPI_type
31 : use m_time, only : timab
32 : use m_pawang, only : pawang_type
33 : use m_pawtab, only : pawtab_type
34 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_allgather, pawcprj_free
35 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle
36 : use m_kg, only : ph1d3d
37 : use m_gsphere, only : symg
38 : use m_fftcore, only : sphereboundary
39 : use m_fft, only : fftpac, fourwf
40 : use m_spacepar, only : symrhg
41 : use m_paw_finegrid, only : pawgylmg
42 : use m_paw_nhat, only : pawsushat
43 :
44 : implicit none
45 :
46 : private
47 : !!***
48 :
49 : public :: suscep_stat ! Compute the susceptibility matrix
50 :
51 : CONTAINS !====================================================================================================
52 : !!***
53 :
54 : !!****f* m_suscep_stat/suscep_stat
55 : !! NAME
56 : !! suscep_stat
57 : !!
58 : !! FUNCTION
59 : !! Compute the susceptibility matrix
60 : !! from input wavefunctions, band occupations, and k point wts.
61 : !! Include the usual sum-over-state terms, but also the
62 : !! corrections due to the change of the Fermi level in the metallic
63 : !! case, as well as implicit sum over higher lying conduction
64 : !! states, thanks to the closure relation (referred to as an extrapolation).
65 : !!
66 : !! INPUTS
67 : !! atindx(natom)=index table for atoms (see scfcv.f)
68 : !! atindx1(natom)=index table for atoms, inverse of atindx (see scfcv.f)
69 : !! cg(2,mcg)=wf in G space
70 : !! cprj(natom,mcg*usecprj)= wave functions projected with non-local projectors:
71 : !! cprj_nk(i)=<p_i|Cnk> where p_i is a non-local projector.
72 : !! dielar(7)=input parameters for dielectric matrix and susceptibility:
73 : !! diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
74 : !! dimcprj(natom*usepaw)=array of dimensions of array cprj (ordered by atom-type)
75 : !! doccde(mband*nkpt*nsppol)=derivative of occupancies wrt
76 : !! the energy for each band and k point
77 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
78 : !! gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for the dielectric matrix
79 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
80 : !! irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
81 : !! istwfk(nkpt)=input option parameter that describes the storage of wfs
82 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
83 : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
84 : !! lmax_diel=1+max. value of l angular momentum used for dielectric matrix
85 : !! mband=maximum number of bands
86 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
87 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
88 : !! mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix
89 : !! mkmem=number of k points treated by this node
90 : !! mpi_enreg=information about MPI parallelization
91 : !! mpw=maximum allowed value for npw
92 : !! natom=number of atoms in cell
93 : !! nband(nkpt*nsppol)=number of bands to be included in summation
94 : !! at each k point for each spin channel
95 : !! neglect_pawhat=1 if PAW contribution from hat density (compensation charge)
96 : !! has to be neglected (to be used when only an estimation of
97 : !! suscep. matrix has to be evaluated, i.e. for SCF precondictioning)
98 : !! nfftdiel=number of fft grid points for the computation of the diel matrix
99 : !! ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix,
100 : !! see ~abinit/doc/variables/vargs.htm#ngfft
101 : !! nkpt=number of k points
102 : !! npwarr(nkpt)=number of planewaves and boundary planewaves
103 : !! at each k, for going from the WF sphere to the medium size FFT grid.
104 : !! npwdiel=third and fifth dimension of the susmat array.
105 : !! nspden=number of spin-density components
106 : !! nspinor=number of spinorial components of the wavefunctions
107 : !! nsppol=1 for unpolarized, 2 for spin-polarized
108 : !! nsym=number of symmetry elements in group (at least 1 for identity)
109 : !! ntypat=number of types of atoms in unit cell.
110 : !! occ(mband*nkpt*nsppol)=
111 : !! occupation numbers for each band (usually 2.0) at each k point
112 : !! occopt=option for occupancies
113 : !! pawang <type(pawang_type)>=paw angular mesh and related data
114 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
115 : !! phnonsdiel(2,nfftdiel**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
116 : !! ph1ddiel(2,3*(2*mgfftdiel+1)*natom*usepaw)=one-dimensional structure factor information
117 : !! for the dielectric matrix
118 : !! rprimd(3,3)=dimensional real space primitive translations
119 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
120 : !! symrel(3,3,nsym)=symmetry matrices in real space (integers)
121 : !! tnons(3,nsym)=reduced nonsymmorphic translations
122 : !! (symrel and tnons are in terms of real space primitive translations)
123 : !! typat(natom)=type (integer) for each atom
124 : !! ucvol=unit cell volume (Bohr**3)
125 : !! unpaw=unit number for cprj PAW data (if used)
126 : !! usecprj= 1 if cprj array is stored in memory
127 : !! usepaw=flag for PAW
128 : !! usetimerev=1 if Time-Reversal symmetry has to be used when symmetrizing susceptibility
129 : !! wtk(nkpt)=k point weights (they sum to 1.0)
130 : !! ylmdiel(npwdiel,lmax_diel**2)= real spherical harmonics for each G and k point
131 : !! for the dielectric matrix
132 : !!
133 : !! OUTPUT
134 : !! susmat(2,npwdiel,nspden,npwdiel,nspden)=
135 : !! the susceptibility (or density-density response) matrix in reciprocal space
136 : !!
137 : !! NOTES
138 : !! Case of non-collinear magnetism:
139 : !! In principle, we should compute 16 susceptibility matrix: chi0-(s1,s2),(s3,s4)
140 : !! (where s1, s2, s3,and s4 are spin indexes)...
141 : !! But, for the time being, the susceptibility is only used to compute the
142 : !! dielectric matrix within RPA approximation; in this approximation, only
143 : !! four susceptibilities are non-zero: chi0-(s1,s1),(s3,s3).
144 : !! They are stored in susmat(:,ipw1,1:2,ipw2,1:2)
145 : !!
146 : !! SOURCE
147 :
148 82 : subroutine suscep_stat(atindx,atindx1,cg,cprj,dielar,dimcprj,doccde,&
149 82 : & eigen,gbound_diel,gprimd,irrzondiel,istwfk,kg,&
150 82 : & kg_diel,lmax_diel,&
151 82 : & mband,mcg,mcprj,mgfftdiel,mkmem,mpi_enreg,mpw,natom,nband,&
152 82 : & neglect_pawhat,nfftdiel,ngfftdiel,nkpt,npwarr,&
153 82 : & npwdiel,nspden,nspinor,nsppol,nsym,ntypat,occ,occopt,&
154 82 : & pawang,pawtab,phnonsdiel,ph1ddiel,rprimd,&
155 82 : & susmat,symafm,symrel,tnons,typat,ucvol,unpaw,usecprj,usepaw,usetimerev,&
156 82 : & wtk,ylmdiel)
157 :
158 : !Arguments ------------------------------------
159 : !scalars
160 : integer,intent(in) :: lmax_diel,mband,mcg,mcprj,mgfftdiel,mkmem,mpw,natom,neglect_pawhat
161 : integer,intent(in) :: nfftdiel,nkpt,npwdiel,nspden,nspinor,nsppol,nsym,ntypat,occopt
162 : integer,intent(in) :: unpaw,usecprj,usepaw,usetimerev
163 : real(dp),intent(in) :: ucvol
164 : type(MPI_type),intent(in) :: mpi_enreg
165 : type(pawang_type),intent(in) :: pawang
166 : !arrays
167 : integer,intent(in) :: atindx(natom),atindx1(natom),dimcprj(natom*usepaw)
168 : integer,intent(in) :: gbound_diel(2*mgfftdiel+8,2)
169 : !no_abirules
170 : !nfftdiel**(1-1/nsym) is 1 if nsym==1, and nfftdiel otherwise
171 : integer,intent(in) :: irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),&
172 : & istwfk(nkpt)
173 : integer,intent(in) :: kg(3,mpw*mkmem),kg_diel(3,npwdiel),&
174 : & nband(nkpt*nsppol),ngfftdiel(18)
175 : integer,intent(in) :: npwarr(nkpt),symafm(nsym),symrel(3,3,nsym),typat(ntypat)
176 : real(dp),intent(in) :: cg(2,mcg),dielar(7)
177 : real(dp),intent(in) :: doccde(mband*nkpt*nsppol),eigen(mband*nkpt*nsppol)
178 : real(dp),intent(in) :: gprimd(3,3),occ(mband*nkpt*nsppol)
179 : !nfftdiel**(1-1/nsym) is 1 if nsym==1, and nfftdiel otherwise
180 : real(dp),intent(in) :: phnonsdiel(2,nfftdiel**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),&
181 : & tnons(3,nsym),wtk(nkpt)
182 : real(dp),intent(in) :: ph1ddiel(2,(3*(2*mgfftdiel+1)*natom)*usepaw),rprimd(3,3)
183 : real(dp),intent(in) :: ylmdiel(npwdiel,lmax_diel**2)
184 : real(dp),intent(out) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
185 : type(pawcprj_type) :: cprj(natom,mcprj*usecprj)
186 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
187 :
188 : !Local variables-------------------------------
189 : !scalars
190 : integer :: bdtot_index,diag,extrap,i1,i2,i3,iband,ibg,icg,ier,ierr
191 : integer :: ifft,ii,ikg,ikpt,indx,iorder_cprj,ipw1,ipw2,isp,isp1,isp2
192 : integer :: ispinor,istwf_k,isym,j1,j2,j3,jj,jsp,k1,k2,k3
193 : integer :: my_nspinor,nband_k,nband_loc,ndiel1,ndiel2,ndiel3,ndiel4,ndiel5,ndiel6
194 : integer :: nkpg_diel,npw_k,npwsp,nspden_eff,nspden_tmp,nsym1,nsym2
195 : integer :: spaceComm,t1,t2,testocc
196 : real(dp) :: ai,ai2,ar,ar2,diegap,dielam,emax,invnsym
197 : real(dp) :: invnsym1,invnsym2,phi1,phi12,phi2,phr1,phr12
198 : real(dp) :: phr2,sumdocc,weight
199 : logical :: antiferro
200 : character(len=500) :: message
201 82 : type(MPI_type) :: mpi_enreg_diel
202 :
203 : !arrays
204 82 : integer,allocatable :: gbound(:,:),kg_k(:,:),sym_g(:,:)
205 82 : integer,allocatable :: tmrev_g(:)
206 : real(dp) :: kpt_diel(3,1),tsec(2)
207 82 : real(dp),allocatable :: drhode(:,:,:),drhode_wk(:,:,:)
208 82 : real(dp),allocatable :: eig_diel(:),gylmg_diel(:,:,:),kpg_dum(:,:)
209 82 : real(dp),allocatable :: occ_deavg(:),ph3d_diel(:,:,:),phdiel(:,:,:)
210 82 : real(dp),allocatable :: phkxred_diel(:,:),rhoextrap(:,:,:,:),rhoextrg(:,:)
211 82 : real(dp),allocatable :: rhoextrr(:,:),sush(:),sussum(:),susvec(:,:,:)
212 82 : real(dp),allocatable :: suswk(:,:,:,:),zhpev1(:,:),zhpev2(:)
213 82 : type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_loc(:,:)
214 :
215 : ! *************************************************************************
216 :
217 82 : call timab(740,1,tsec)
218 82 : call timab(741,1,tsec)
219 :
220 82 : ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
221 :
222 :
223 : !----- Initialisations -----------------------------------------------------------
224 : !---------------------------------------------------------------------------------
225 :
226 82 : if (usecprj==0.and.usepaw==1) then
227 : write (message,'(3a)')&
228 0 : & ' cprj datastructure must be allocated !',ch10,&
229 0 : & ' Action: change pawusecp input keyword.'
230 0 : ABI_ERROR(message)
231 : end if
232 :
233 82 : if (mpi_enreg%paral_spinor==1) then
234 0 : message = ' not yet allowed for parallelization over spinors !'
235 0 : ABI_ERROR(message)
236 : end if
237 :
238 : !Init mpicomm
239 82 : if(mpi_enreg%paral_kgb==1) then
240 0 : spaceComm=mpi_enreg%comm_kpt
241 : else
242 82 : spaceComm=mpi_enreg%comm_cell
243 : end if
244 :
245 : !The dielectric stuff is performed in sequential mode.
246 : !Set mpi_enreg_diel accordingly
247 82 : call initmpi_seq(mpi_enreg_diel)
248 82 : call MPI_enreg_diel%distribfft%init_seq('c',ngfftdiel(2),ngfftdiel(3),'all')
249 :
250 : !testocc to be taken away
251 82 : testocc=1
252 :
253 82 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
254 :
255 82 : iorder_cprj=0 ! order for the cprj reading...
256 :
257 : !Initialize some scalar quantities
258 82 : antiferro=(nsppol==1.and.nspden==2)
259 82 : nspden_eff=min(max(nsppol,nspden),2) ! Size for the computed part of susmat
260 82 : bdtot_index=0 ; icg=0 ; ibg=0
261 82 : ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3)
262 :
263 : !ndiel4,ndiel5,ndiel6 are FFT dimensions, modified to avoid cache trashing
264 82 : ndiel4=ngfftdiel(4) ; ndiel5=ngfftdiel(5) ; ndiel6=ngfftdiel(6)
265 82 : diegap=dielar(5) ; dielam=dielar(6)
266 82 : extrap=0
267 :
268 : !If dielam is too small, there is no extrapolation.
269 82 : if(dielam>1.0d-6)extrap=1
270 :
271 : !Some stuff for symmetries
272 1006 : nsym1=sum(symafm,mask=symafm==1)
273 82 : nsym2=nsym-nsym1
274 82 : invnsym =one/dble(nsym)
275 82 : invnsym1=one/dble(nsym1)
276 82 : invnsym2=one
277 : !FIXME: make sure this is consistent with following code
278 : !div by 0 for several v5 tests
279 82 : if (nsym2 > 0) invnsym2=one/dble(nsym2)
280 :
281 : !Allocations
282 246 : ABI_MALLOC(occ_deavg,(mband))
283 82 : if(occopt>=3) then
284 228 : ABI_MALLOC(drhode,(2,npwdiel,nspden_eff))
285 : else
286 25 : ABI_MALLOC(drhode,(0,0,0))
287 : end if
288 82 : if(extrap==1) then
289 486 : ABI_MALLOC(rhoextrap,(ndiel4,ndiel5,ndiel6,nspinor))
290 : else
291 82 : ABI_MALLOC(rhoextrap,(0,0,0,0))
292 : end if
293 :
294 : !zero the susceptibility matrix and other needed quantities
295 8744382 : susmat(:,:,:,:,:)=zero
296 82 : if(occopt>=3)then
297 16745 : drhode(:,:,:)=zero
298 57 : sumdocc=zero
299 : end if
300 :
301 : !PAW additional initializations
302 82 : if (usepaw==1) then
303 0 : ABI_MALLOC(gylmg_diel,(npwdiel,lmax_diel**2,ntypat))
304 0 : ABI_MALLOC(ph3d_diel,(2,npwdiel,natom))
305 0 : if (neglect_pawhat==0) then
306 0 : ABI_MALLOC(phkxred_diel,(2,natom))
307 0 : ABI_MALLOC(kpg_dum,(0,0))
308 0 : kpt_diel(1:3,1)=zero;phkxred_diel(1,:)=one;phkxred_diel(2,:)=zero;nkpg_diel=0
309 : ! write(std_out,*) ' lmax_diel ', lmax_diel
310 0 : call pawgylmg(gprimd,gylmg_diel,kg_diel,kpg_dum,kpt_diel,lmax_diel,nkpg_diel,npwdiel,ntypat,pawtab,ylmdiel)
311 0 : call ph1d3d(1,natom,kg_diel,natom,natom,npwdiel,ndiel1,ndiel2,ndiel3,phkxred_diel,ph1ddiel,ph3d_diel)
312 0 : ABI_FREE(phkxred_diel)
313 0 : ABI_FREE(kpg_dum)
314 : else
315 0 : gylmg_diel=zero;ph3d_diel=one
316 : end if
317 : else
318 82 : ABI_MALLOC(gylmg_diel,(0,0,0))
319 82 : ABI_MALLOC(ph3d_diel,(0,0,0))
320 : end if
321 :
322 82 : call timab(741,2,tsec)
323 :
324 :
325 :
326 : !--BIG loop over spins ------------------------------------------------------------
327 : !---------------------------------------------------------------------------------
328 :
329 168 : do isp=1,nsppol
330 86 : ikg=0
331 :
332 307385 : if(extrap==1)rhoextrap(:,:,:,:)=zero
333 :
334 : ! --BIG loop over k-points --------------------------------------------------------
335 : ! ---------------------------------------------------------------------------------
336 :
337 424 : do ikpt=1,nkpt
338 :
339 338 : nband_k=nband(ikpt+(isp-1)*nkpt)
340 338 : istwf_k=istwfk(ikpt)
341 338 : npw_k=npwarr(ikpt)
342 :
343 338 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isp,mpi_enreg%me_kpt)) then
344 16 : bdtot_index=bdtot_index+nband_k
345 16 : cycle
346 : end if
347 :
348 322 : call timab(742,1,tsec)
349 :
350 :
351 1610 : ABI_MALLOC(gbound,(2*mgfftdiel+8,2))
352 966 : ABI_MALLOC(kg_k,(3,npw_k))
353 :
354 322 : if (usepaw==1) then
355 0 : ABI_MALLOC(cprj_k,(natom,my_nspinor*nband_k))
356 0 : if (neglect_pawhat==0) then
357 0 : call pawcprj_alloc(cprj_k,0,dimcprj)
358 0 : if (mpi_enreg%nproc_band==1) then
359 : call pawcprj_get(atindx1,cprj_k,cprj,natom,1,ibg,ikpt,iorder_cprj,isp,&
360 : & mband,mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,unpaw,&
361 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
362 : else
363 0 : nband_loc=nband_k/mpi_enreg%nproc_band
364 0 : ABI_MALLOC(cprj_loc,(natom,my_nspinor*nband_loc))
365 0 : call pawcprj_alloc(cprj_loc,0,dimcprj)
366 : call pawcprj_get(atindx1,cprj_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isp,&
367 : & mband/mpi_enreg%nproc_band,mkmem,natom,nband_loc,nband_loc,my_nspinor,nsppol,unpaw,&
368 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
369 : call pawcprj_mpi_allgather(cprj_loc,cprj_k,natom,my_nspinor*nband_loc,mpi_enreg%bandpp,&
370 0 : & dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.true.)
371 0 : call pawcprj_free(cprj_loc)
372 0 : ABI_FREE(cprj_loc)
373 : end if
374 : else
375 : !call pawcprj_nullify(cprj_k)
376 : end if
377 : else
378 322 : ABI_MALLOC(cprj_k,(0,0))
379 : end if
380 :
381 562306 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
382 322 : call sphereboundary(gbound,istwf_k,kg_k,mgfftdiel,npw_k)
383 :
384 322 : if(extrap==1)then
385 : ! Compute inverse of average dielectric gap for each band
386 : ! and multiply by occupation factor
387 4598 : emax=maxval(eigen(1+bdtot_index:nband_k+bdtot_index))
388 4277 : do iband=1,nband_k
389 : occ_deavg(iband)= occ(iband+bdtot_index)*dielam &
390 4277 : & / ( emax-eigen(iband+bdtot_index) + diegap )
391 : end do
392 : else
393 21 : occ_deavg(:)=zero
394 : end if
395 :
396 322 : call timab(742,2,tsec)
397 322 : call timab(743,1,tsec)
398 :
399 : ! Compute the contribution of each k-point to susmat, rhoextrap, drhode and sumdocc.
400 322 : if(mpi_enreg%paral_kgb==1)then !Only this version is in parallel
401 : ! Use either the simpler implementation
402 : ! Should provide a test !!!
403 : call susk(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,&
404 : & gbound_diel,gylmg_diel,icg,ikpt,&
405 : & isp,istwf_k,kg_diel,kg_k,lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,&
406 : & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,&
407 : & npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,&
408 : & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,&
409 0 : & susmat,typat,ucvol,usepaw,wtk)
410 : else
411 : ! Or the more sophisticated one, needed to save memory.
412 : call suskmm(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,&
413 : & gbound_diel,gylmg_diel,icg,ikpt,&
414 : & isp,istwf_k,kg_diel,kg_k,lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,&
415 : & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,&
416 : & npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,&
417 : & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,&
418 322 : & susmat,typat,ucvol,usepaw,wtk)
419 : end if
420 :
421 322 : call timab(743,2,tsec)
422 :
423 322 : ABI_FREE(gbound)
424 322 : ABI_FREE(kg_k)
425 :
426 322 : bdtot_index=bdtot_index+nband_k
427 :
428 322 : if (mkmem/=0) then
429 322 : ibg=ibg+my_nspinor*nband_k
430 322 : icg=icg+my_nspinor*npw_k*nband_k
431 322 : ikg=ikg+npw_k
432 : end if
433 322 : if (usepaw==1) then
434 0 : if (neglect_pawhat==0) then
435 0 : call pawcprj_free(cprj_k)
436 : end if
437 : end if
438 408 : ABI_FREE(cprj_k)
439 :
440 : ! End loop on ikpt: --------------------------------------------------------
441 : end do
442 :
443 : ! Here include the contribution from the extrapolation to susmat,
444 : ! diagonal part
445 168 : if(extrap==1)then
446 :
447 85 : call timab(744,1,tsec)
448 :
449 : ! Transfer extrapolating density on augmented fft grid to
450 : ! normal fft grid in real space.
451 : ! Warning1 : if collinear magnetism, must treat only one spin at a time
452 : ! Warning2 : if non-collinear magnetism, treat both spins
453 : ! Warning3 : this is subtle for antiferro magnetism
454 85 : nspden_tmp=1;if (antiferro) nspden_tmp=2
455 340 : ABI_MALLOC(rhoextrr,(nfftdiel,nspden_tmp))
456 255 : ABI_MALLOC(rhoextrg,(2,nfftdiel))
457 85 : if (nspden==1.and.nspinor==2) rhoextrap(:,:,:,1)=rhoextrap(:,:,:,1)+rhoextrap(:,:,:,2)
458 :
459 170 : do ispinor=1,min(nspinor,nspden)
460 85 : jsp=isp+ispinor-1
461 :
462 : call fftpac(1,mpi_enreg_diel,1,ndiel1,ndiel2,ndiel3,ndiel4,ndiel5,ndiel6,&
463 85 : & ngfftdiel,rhoextrr(:,1),rhoextrap(:,:,:,ispinor),1)
464 :
465 : ! Generate the density in reciprocal space, and symmetrize it
466 : ! (note symrhg also make the reverse FFT, to get symmetrized density;
467 : ! this is useless here, and should be made an option)
468 : call symrhg(1,gprimd,irrzondiel,mpi_enreg_diel,nfftdiel,nfftdiel,ngfftdiel,&
469 85 : & nspden_tmp,1,nsym,phnonsdiel,rhoextrg,rhoextrr,rprimd,symafm,symrel,tnons)
470 :
471 12559 : do ipw2=1,npwdiel
472 12389 : j1=kg_diel(1,ipw2) ; j2=kg_diel(2,ipw2) ; j3=kg_diel(3,ipw2)
473 : ! static: Only fills lower half of the matrix (here, the susceptibility matrix)
474 : ! dynamical: fill all, will not affect susopt==2 for which extrap==0
475 2872583 : do ipw1=1,npwdiel
476 2860109 : i1=kg_diel(1,ipw1) ; i2=kg_diel(2,ipw1) ; i3=kg_diel(3,ipw1)
477 : ! NOTE that there is a FFT folding (superposition) bias here
478 : ! Should use kgindex, in the same spirit as in prcref
479 2860109 : k1=i1-j1; k1=modulo(k1,ndiel1)
480 2860109 : k2=i2-j2; k2=modulo(k2,ndiel2)
481 2860109 : k3=i3-j3; k3=modulo(k3,ndiel3)
482 2860109 : ifft=k1+1+ndiel1*(k2+ndiel2*k3)
483 2860109 : susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+rhoextrg(1,ifft)
484 2872498 : susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+rhoextrg(2,ifft)
485 : end do
486 : end do
487 :
488 : end do
489 85 : ABI_FREE(rhoextrg)
490 85 : ABI_FREE(rhoextrr)
491 :
492 85 : call timab(744,2,tsec)
493 :
494 : end if
495 :
496 : ! End loop over spins ---------------------------------------------------------
497 : end do
498 :
499 82 : ABI_FREE(occ_deavg)
500 82 : ABI_FREE(rhoextrap)
501 82 : ABI_FREE(gylmg_diel)
502 82 : ABI_FREE(ph3d_diel)
503 : !end if
504 :
505 82 : call destroy_mpi_enreg(mpi_enreg_diel)
506 :
507 : !-- Stuff for parallelism --------------------------------------------------------
508 : !---------------------------------------------------------------------------------
509 :
510 : if(xmpi_paral==1)then
511 82 : call timab(746,1,tsec)
512 246 : ABI_MALLOC(sussum,(2*npwdiel*nspden*npwdiel*nspden))
513 : ! Recreate full susmat on all proc.
514 : ! This should be coded more efficiently,
515 : ! since half of the matrix is still empty, and
516 : ! it is spin-diagonal.
517 164 : sussum(:)=reshape(susmat(:,:,:,:,:),(/2*npwdiel*nspden*npwdiel*nspden/))
518 82 : call xmpi_sum(sussum,spaceComm,ierr)
519 492 : susmat(:,:,:,:,:)=reshape(sussum(:),(/2,npwdiel,nspden,npwdiel,nspden/))
520 82 : ABI_FREE(sussum)
521 : ! Recreate full drhode on all proc.
522 82 : if(occopt>=3 .and. testocc==1)then
523 57 : call xmpi_sum(drhode,spaceComm,ierr)
524 : ! Should use only one mpi-allreduce call instead of the three
525 57 : call xmpi_sum(sumdocc,spaceComm,ierr)
526 : end if
527 164 : call timab(746,2,tsec)
528 : end if
529 :
530 : !-- Apply spatial hermitian/symmetries on spin-diagonal susceptibility matrix ----
531 : !---------------------------------------------------------------------------------
532 :
533 82 : call timab(747,1,tsec)
534 :
535 : !If antiferro magnetism, has to divide (spin-diagonal) susceptibility by 2 (due to dble occupations)
536 82 : if (antiferro) then
537 0 : do ipw2=1,npwdiel
538 0 : do ipw1=ipw2,npwdiel
539 0 : susmat(:,ipw1,1,ipw2,1)=half*susmat(:,ipw1,1,ipw2,1)
540 : end do
541 : end do
542 : end if
543 :
544 : !Generate upper half of the spin-diagonal matrix (still the susceptibility matrix)
545 168 : do isp=1,nspden_eff
546 12648 : do ipw2=2,npwdiel
547 1452002 : do ipw1=1,ipw2-1
548 1439436 : susmat(1,ipw1,isp,ipw2,isp)= susmat(1,ipw2,isp,ipw1,isp)
549 1451916 : susmat(2,ipw1,isp,ipw2,isp)=-susmat(2,ipw2,isp,ipw1,isp)
550 : end do
551 : end do
552 : end do
553 :
554 : !Compute symmetric of G-vectors and eventual phases
555 : !(either time-reversal or spatial symmetries)
556 246 : ABI_MALLOC(tmrev_g,(npwdiel))
557 328 : ABI_MALLOC(sym_g,(npwdiel,nsym))
558 328 : ABI_MALLOC(phdiel,(2,npwdiel,nsym))
559 82 : call symg(kg_diel,npwdiel,nsym,phdiel,sym_g,symrel,tmrev_g,tnons)
560 :
561 : !Impose spatial symmetries to the spin-diagonal susceptibility matrix
562 410 : ABI_MALLOC(suswk,(2,npwdiel,npwdiel,nspden_eff))
563 168 : do isp=1,nspden_eff
564 8687048 : suswk(:,:,:,isp)=susmat(:,:,isp,:,isp) ! Temporary storage
565 : end do
566 :
567 168 : do isp=1,nspden_eff
568 86 : jsp=min(3-isp,nsppol)
569 12734 : do ipw2=1,npwdiel
570 2904090 : do ipw1=1,npwdiel
571 2891438 : ar=suswk(1,ipw1,ipw2,isp)
572 2891438 : ai=suswk(2,ipw1,ipw2,isp)
573 2891438 : ar2=zero;ai2=zero
574 2891438 : if(nsym>1)then
575 8888932 : do isym=2,nsym
576 8230736 : t1=sym_g(ipw1,isym) ; t2=sym_g(ipw2,isym)
577 : ! Not all symmetries are non-symmorphic. Should save time here ...
578 8230736 : phr1=phdiel(1,ipw1,isym) ; phi1=phdiel(2,ipw1,isym)
579 8230736 : phr2=phdiel(1,ipw2,isym) ; phi2=phdiel(2,ipw2,isym)
580 8230736 : phr12= phr1*phr2+phi1*phi2 ; phi12=phi1*phr2-phr1*phi2
581 8888932 : if (symafm(isym)==1) then
582 8230736 : ar=ar+suswk(1,t1,t2,isp)*phr12-suswk(2,t1,t2,isp)*phi12
583 8230736 : ai=ai+suswk(2,t1,t2,isp)*phr12+suswk(1,t1,t2,isp)*phi12
584 : else
585 0 : ar2=ar2+suswk(1,t1,t2,jsp)*phr12-suswk(2,t1,t2,jsp)*phi12
586 0 : ai2=ai2+suswk(2,t1,t2,jsp)*phr12+suswk(1,t1,t2,jsp)*phi12
587 : end if
588 : end do
589 : end if
590 2904004 : if (antiferro) then
591 0 : susmat(1,ipw1,1,ipw2,1)=ar*invnsym1
592 0 : susmat(2,ipw1,1,ipw2,1)=ai*invnsym1
593 0 : susmat(1,ipw1,2,ipw2,2)=ar2*invnsym2
594 0 : susmat(2,ipw1,2,ipw2,2)=ai2*invnsym2
595 : else
596 2891438 : susmat(1,ipw1,isp,ipw2,isp)=(ar+ar2)*invnsym
597 2891438 : susmat(2,ipw1,isp,ipw2,isp)=(ai+ai2)*invnsym
598 : end if
599 : end do
600 : end do
601 : end do
602 82 : ABI_FREE(suswk)
603 :
604 :
605 : !-- Add contribibution to susceptibility due to change of Fermi level -----------
606 : !---------------------------------------------------------------------------------
607 :
608 82 : if (occopt>=3.and.testocc==1) then
609 :
610 : ! Impose spatial symmetries to drhode
611 228 : ABI_MALLOC(drhode_wk,(2,npwdiel,nspden_eff))
612 116 : do isp=1,nspden_eff
613 59 : jsp=min(3-isp,nsppol)
614 5659 : do ipw1=1,npwdiel
615 5543 : ar=drhode(1,ipw1,isp)
616 5543 : ai=drhode(2,ipw1,isp)
617 5543 : ar2=zero;ai2=zero
618 5543 : if (nsym>1) then
619 79744 : do isym=2,nsym
620 74201 : t1=sym_g(ipw1,isym)
621 : ! Not all symmetries are non-symmorphic. Should save time here ...
622 74201 : phr1=phdiel(1,ipw1,isym);phi1=phdiel(2,ipw1,isym)
623 79744 : if (symafm(isym)==1) then
624 74201 : ar=ar+drhode(1,t1,isp)*phr1-drhode(2,t1,isp)*phi1
625 74201 : ai=ai+drhode(2,t1,isp)*phr1+drhode(1,t1,isp)*phi1
626 : else
627 0 : ar2=ar2+drhode(1,t1,jsp)*phr1-drhode(2,t1,jsp)*phi1
628 0 : ai2=ai2+drhode(2,t1,jsp)*phr1+drhode(1,t1,jsp)*phi1
629 : end if
630 : end do
631 : end if
632 5602 : if (antiferro) then ! 1/2 factor due to (dble) occupations
633 0 : drhode_wk(1,ipw1,1)=half*ar*invnsym1
634 0 : drhode_wk(2,ipw1,1)=half*ai*invnsym1
635 0 : drhode_wk(1,ipw1,2)=half*ar2*invnsym2
636 0 : drhode_wk(2,ipw1,2)=half*ai2*invnsym2
637 : else
638 5543 : drhode_wk(1,ipw1,isp)=(ar+ar2)*invnsym
639 5543 : drhode_wk(2,ipw1,isp)=(ai+ai2)*invnsym
640 : end if
641 : end do
642 : end do
643 :
644 : ! Add contribution to non-diagonal susceptibility
645 : ! Presently fills complete susceptibility matrix, not only lower half
646 57 : weight=one/sumdocc
647 116 : do isp2=1,nspden_eff
648 5659 : do ipw2=1,npwdiel
649 11205 : do isp1=1,nspden_eff
650 638913 : do ipw1=1,npwdiel
651 : susmat(1,ipw1,isp1,ipw2,isp2)=susmat(1,ipw1,isp1,ipw2,isp2)- &
652 : & weight*( drhode_wk(1,ipw1,isp1)*drhode_wk(1,ipw2,isp2) &
653 627767 : & +drhode_wk(2,ipw1,isp1)*drhode_wk(2,ipw2,isp2) )
654 : susmat(2,ipw1,isp1,ipw2,isp2)=susmat(2,ipw1,isp1,ipw2,isp2)- &
655 : & weight*( drhode_wk(2,ipw1,isp1)*drhode_wk(1,ipw2,isp2) &
656 633370 : & -drhode_wk(1,ipw1,isp1)*drhode_wk(2,ipw2,isp2) )
657 : end do
658 : end do
659 : end do
660 : end do
661 57 : ABI_FREE(drhode_wk)
662 :
663 : end if
664 : !if (occopt>=3) then
665 82 : ABI_FREE(drhode)
666 : !end if
667 :
668 :
669 : !--- Impose the time-reversal symmetry to the susceptibility matrix --------------
670 : !---------------------------------------------------------------------------------
671 :
672 82 : if (usetimerev==1) then
673 328 : ABI_MALLOC(suswk,(2,npwdiel,npwdiel,1))
674 :
675 : ! Impose the time-reversal symmetry to the spin-diagonal susceptibility matrix
676 168 : do isp=1,nspden_eff
677 8686966 : suswk(:,:,:,1)=susmat(:,:,isp,:,isp) ! Temporary storage
678 12734 : do ipw2=1,npwdiel
679 12566 : t2=tmrev_g(ipw2)
680 2904090 : do ipw1=1,npwdiel
681 2891438 : t1=tmrev_g(ipw1)
682 2891438 : susmat(1,ipw1,isp,ipw2,isp)=half*(suswk(1,ipw1,ipw2,1)+suswk(1,t1,t2,1))
683 2904004 : susmat(2,ipw1,isp,ipw2,isp)=half*(suswk(2,ipw1,ipw2,1)-suswk(2,t1,t2,1))
684 : end do
685 : end do
686 : end do
687 :
688 : ! Impose the time-reversal symmetry to the off-diagonal susceptibility matrix
689 82 : if (nspden_eff/=1.and.occopt>=3.and.testocc==1) then
690 1382 : suswk(:,:,:,1)=susmat(:,:,1,:,2) ! Temporary storage
691 32 : do ipw2=1,npwdiel
692 30 : t2=tmrev_g(ipw2)
693 482 : do ipw1=1,npwdiel
694 450 : t1=tmrev_g(ipw1)
695 450 : ar=half*(suswk(1,ipw1,ipw2,1)+suswk(1,t1,t2,1))
696 450 : ai=half*(suswk(2,ipw1,ipw2,1)-suswk(2,t1,t2,1))
697 450 : susmat(1,ipw1,1,ipw2,2)= ar
698 450 : susmat(2,ipw1,1,ipw2,2)= ai
699 450 : susmat(1,ipw1,2,ipw2,1)= ar
700 480 : susmat(2,ipw1,2,ipw2,1)=-ai
701 : end do
702 : end do
703 : end if
704 82 : ABI_FREE(suswk)
705 : end if
706 :
707 82 : ABI_FREE(phdiel)
708 82 : ABI_FREE(sym_g)
709 82 : ABI_FREE(tmrev_g)
710 :
711 :
712 : !-- The full susceptibility matrix is computed -----------------------------------
713 : !-- Now, eventually diagonalize it and stop --------------------------------------
714 : !---------------------------------------------------------------------------------
715 :
716 : !Must turn on this flag to make the diagonalisation
717 82 : diag=0
718 : if(diag==1)then
719 :
720 : npwsp=npwdiel*nspden_eff
721 : ABI_MALLOC(sush,(npwsp*(npwsp+1)))
722 : ABI_MALLOC(susvec,(2,npwsp,npwsp))
723 : ABI_MALLOC(eig_diel,(npwsp))
724 : ABI_MALLOC(zhpev1,(2,2*npwsp-1))
725 : ABI_MALLOC(zhpev2,(3*npwsp-2))
726 : ier=0
727 :
728 : ! Store the susceptibility matrix in proper mode before calling zhpev
729 : indx=1
730 : do ii=1,npwdiel
731 : do jj=1,ii
732 : sush(indx )=susmat(1,jj,1,ii,1)
733 : sush(indx+1)=susmat(2,jj,1,ii,1)
734 : indx=indx+2
735 : end do
736 : end do
737 :
738 : ! If spin-polarized, need to store other parts of the matrix
739 : if(nspden_eff/=1)then
740 : do ii=1,npwdiel
741 : ! Here, spin-flip contribution
742 : do jj=1,npwdiel
743 : sush(indx )=susmat(1,jj,1,ii,2)
744 : sush(indx+1)=susmat(2,jj,1,ii,2)
745 : indx=indx+2
746 : end do
747 : ! Here spin down-spin down upper matrix
748 : do jj=1,ii
749 : sush(indx )=susmat(1,jj,2,ii,2)
750 : sush(indx+1)=susmat(2,jj,2,ii,2)
751 : indx=indx+2
752 : end do
753 : end do
754 : end if
755 :
756 : call ZHPEV ('V','U',npwsp,sush,eig_diel,susvec,npwsp,zhpev1,&
757 : & zhpev2,ier)
758 :
759 : write(std_out,*)' suscep_stat : print eigenvalues of the susceptibility matrix'
760 : do ii=1,npwsp
761 : write(std_out,'(i5,es16.6)' )ii,eig_diel(ii)
762 : end do
763 :
764 : ABI_FREE(sush)
765 : ABI_FREE(susvec)
766 : ABI_FREE(eig_diel)
767 : ABI_FREE(zhpev1)
768 : ABI_FREE(zhpev2)
769 : ABI_ERROR("Stopping here!")
770 : end if
771 :
772 82 : call timab(747,2,tsec)
773 82 : call timab(740,2,tsec)
774 :
775 82 : end subroutine suscep_stat
776 : !!***
777 :
778 : !!****f* m_suscep_stat/susk
779 : !! NAME
780 : !! susk
781 : !!
782 : !! FUNCTION
783 : !! Compute the contribution of one k point to the susceptibility matrix
784 : !! from input wavefunctions, band occupations, and k point wts.
785 : !! Include the usual sum-over-state terms, but also the
786 : !! corrections due to the change of the Fermi level in the metallic
787 : !! case, as well as implicit sum over higher lying conduction
788 : !! states, thanks to the closure relation (referred to as an extrapolation).
789 : !! Compared to the routine suskmm, there is no particular attention
790 : !! to the use of the memory, so the code is simpler.
791 : !!
792 : !! INPUTS
793 : !! atindx(natom)=index table for atoms
794 : !! bdtot_index=index for the number of the band
795 : !! cg(2,mcg)=wfs in G space
796 : !! cprj_k(natom,nspinor*nband_k)= wave functions projected with non-local projectors:
797 : !! cprj_k=<p_i|Cnk> where p_i is a non-local projector.
798 : !! doccde(mband*nkpt*nsppol)=derivative of occupancies wrt
799 : !! the energy for each band and k point
800 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
801 : !! extrap: if==1, the closure relation (an extrapolation) must be used
802 : !! gbound(2*mgfftdiel+8,2)=G sphere boundary for going from WF sphere to
803 : !! medium size FFT grid
804 : !! gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for going from medium size
805 : !! FFT grid to small sphere.
806 : !! gylmg_diel(npwdiel,lmax_diel,ntypat*usepaw)= -PAW only- Fourier transform of g_l(r).Y_ml(r) shape functions
807 : !! for dielectric matrix
808 : !! icg=index for cg
809 : !! ikpt=number of the k point
810 : !! isp=number of the current spin
811 : !! istwf_k=input option parameter that describes the storage of wfs
812 : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
813 : !! kg_k(3,npw_k)=coordinates of planewaves in basis sphere.
814 : !! lmax_diel=1+max. value of l angular momentum used for dielectric matrix
815 : !! mband=maximum number of bands
816 : !! mcg=dimension of cg
817 : !! mgfftdiel=maximum size of 1D FFTs, for the computation of
818 : !! the dielectric matrix
819 : !! mpi_enreg=information about MPI parallelization
820 : !! natom=number of atoms in cell
821 : !! nband_k=number of bands at this k point for that spin polarization
822 : !! ndiel4,ndiel5,ndiel6= FFT dimensions, modified to avoid cache trashing
823 : !! neglect_pawhat=1 if PAW contribution from hat density (compensation charge)
824 : !! has to be neglected (to be used when only an estimation of
825 : !! suscep. matrix has to be evaluated, i.e. for SCF precondictioning)
826 : !! ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix,
827 : !! see ~abinit/doc/variables/vargs.htm#ngfft
828 : !! nkpt=number of k points
829 : !! npwdiel=third and fifth dimension of the susmat array.
830 : !! npw_k=number of plane waves at this k point
831 : !! nspden=number of spin-density components
832 : !! nspden_eff=number of spin-density components actually computed in sussceptibility
833 : !! nspinor=number of spinorial components of the wavefunctions
834 : !! nsppol=1 for unpolarized, 2 for spin-polarized
835 : !! ntypat=number of types of atoms in unit cell.
836 : !! occ(mband*nkpt*nsppol)=
837 : !! occupation numbers for each band (usually 2.0) at each k point
838 : !! occopt=option for occupancies
839 : !! occ_deavg(mband)=factor for extrapolation (occup. divided by an energy gap)
840 : !! pawang <type(pawang_type)>=paw angular mesh and related data
841 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
842 : !! ph3d_diel(2,npwdiel,natom*usepaw)=3-dim structure factors, for each atom and plane wave, for dielectric matrix
843 : !! typat(natom)=type (integer) for each atom
844 : !! ucvol=unit cell volume (Bohr**3)
845 : !! usepaw=flag for PAW
846 : !! wtk(nkpt)=k point weights (they sum to 1.0)
847 : !!
848 : !! OUTPUT
849 : !! (see side effects)
850 : !!
851 : !! SIDE EFFECTS
852 : !! These quantities are accumulated in this routine:
853 : !! drhode(2,npwdiel,nspden_eff)=weighted density, needed to compute the
854 : !! effect of change of fermi energy
855 : !! rhoextrap(ndiel4,ndiel5,ndiel6,nspinor)=density-like array, needed for the
856 : !! extrapolation procedure.
857 : !! sumdocc=sum of weighted occupation numbers, needed to compute the
858 : !! effect of change of fermi energy
859 : !! susmat(2,npwdiel,nspden,npwdiel,nspden)=
860 : !! the susceptibility (or density-density response) matrix in reciprocal space
861 : !!
862 : !! NOTES
863 : !! Band-fft parallel treatment: Each processor will treat his own band, but susmat will be known by all.
864 : !! This means that cg will not have the same meaning in sequential or parallel mode.
865 : !! In parallel mode, it will contain the set of all bands treated by the currrent processor.
866 : !! To achieve this, the argument cg has been replaced by cg_mpi, with the "target" attribute.
867 : !! In sequential mode, the pointer cg will point towards cg_mpi. In parallel mode, cg will point
868 : !! to a new array cg_local, containing the bands treated by the currrent processor.
869 : !! This allows to minimize the overhead incurred by the parallelization of the sequential version.
870 : !! A similar treatment is performed on kg_k, npw_k.
871 : !! A future version might have objects like kg_k_gather as arguments, instead of computing them.
872 : !! This is in slight violation of programming rules, but I think it is safe, since the pointers remain local
873 : !! GZ
874 : !!
875 : !! SOURCE
876 :
877 0 : subroutine susk(atindx,bdtot_index,cg_mpi,cprj_k,doccde,drhode,eigen,extrap,gbound,&
878 0 : & gbound_diel,gylmg_diel,icg_mpi,ikpt,isp,istwf_k,kg_diel,kg_k_mpi,&
879 : & lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,&
880 : & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,&
881 0 : & npwdiel,npw_k_mpi,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,&
882 0 : & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,&
883 0 : & susmat,typat,ucvol,usepaw,wtk)
884 :
885 : !Arguments ------------------------------------
886 : !This type is defined in defs_mpi
887 : !scalars
888 : integer,intent(in) :: bdtot_index,extrap,ikpt,isp,istwf_k,lmax_diel,mband,mcg
889 : integer,intent(in) :: mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat
890 : integer,intent(in) :: nkpt,npwdiel,nspden,nspden_eff,nspinor,nsppol
891 : integer,intent(in) :: ntypat,occopt,usepaw
892 : integer,intent(in),target :: icg_mpi,npw_k_mpi
893 : real(dp),intent(in) :: ucvol
894 : real(dp),intent(inout) :: sumdocc
895 : type(MPI_type),intent(in) :: mpi_enreg
896 : type(pawang_type),intent(in) :: pawang
897 : !arrays
898 : integer,intent(in) :: atindx(natom),gbound_diel(2*mgfftdiel+8,2)
899 : integer,intent(in) :: kg_diel(3,npwdiel),ngfftdiel(18),typat(natom)
900 : integer,intent(in),target :: kg_k_mpi(3,npw_k_mpi)
901 : integer,intent(inout) :: gbound(2*mgfftdiel+8,2)
902 0 : integer,pointer :: kg_k(:,:)
903 : real(dp),intent(in) :: doccde(mband*nkpt*nsppol),eigen(mband*nkpt*nsppol)
904 : real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat*usepaw)
905 : real(dp),intent(in) :: occ(mband*nkpt*nsppol),occ_deavg(mband)
906 : real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom*usepaw),wtk(nkpt)
907 : real(dp),intent(in),target :: cg_mpi(2,mcg)
908 : real(dp),intent(inout) :: drhode(2,npwdiel,nspden_eff)
909 : real(dp),intent(inout) :: rhoextrap(ndiel4,ndiel5,ndiel6,nspinor)
910 : real(dp),intent(inout) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
911 : type(pawcprj_type) :: cprj_k(natom,nspinor*nband_k*usepaw)
912 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
913 :
914 : !Local variables-------------------------------
915 : ! real(dp), allocatable :: cg_disk(:,:)
916 : !Local variables for MPI
917 : !scalars
918 : integer :: blocksize,i1,i2,i3,iband,iband_loc,ibd1,ibd2,ibdblock,ier
919 : integer :: iproc,iproc_fft,ipw,ipw1,ipw2,isp1,isp2,ispinor,iwf,jsp,me_bandfft
920 : integer :: nbdblock,ndatarecv,ndiel1,ndiel2,ndiel3
921 : integer :: sizemax_per_proc,spaceComm,testocc,tim_fourwf
922 : integer,pointer :: icg,npw_k
923 : integer,target :: icg_loc=0,npw_k_loc,npw_tot
924 : real(dp) :: eigdiff,occdiff,tolocc,weight,wght1,wght2
925 0 : type(MPI_type) :: mpi_enreg_diel
926 : !arrays
927 0 : integer,allocatable :: band_loc(:),kg_k_gather(:,:),npw_per_proc(:),rdispls(:)
928 0 : integer,allocatable :: rdispls_all(:),rdisplsloc(:),recvcounts(:)
929 0 : integer,allocatable :: recvcountsloc(:),sdispls(:),sdisplsloc(:),sendcounts(:)
930 0 : integer,allocatable :: sendcountsloc(:)
931 0 : integer,allocatable,target :: kg_k_gather_all(:,:)
932 : real(dp) :: tsec(2)
933 0 : real(dp),allocatable :: cwavef(:,:),cwavef_alltoall(:,:)
934 0 : real(dp),allocatable :: cwavef_alltoall_gather(:,:),dummy(:,:),rhoaug(:,:,:)
935 0 : real(dp),allocatable :: susmat_mpi(:,:,:)
936 0 : real(dp),allocatable :: wfprod(:,:),wfraug(:,:,:,:),wfrspa(:,:,:,:,:,:)
937 0 : real(dp),allocatable,target :: cg_local(:,:)
938 0 : real(dp),pointer :: cg(:,:)
939 0 : logical,allocatable :: treat_band(:)
940 :
941 : ! *************************************************************************
942 :
943 : !DEBUG
944 : !write(std_out,*)' susk : enter '; stop
945 : !ENDDEBUG
946 :
947 0 : call timab(750,1,tsec)
948 0 : call timab(751,1,tsec)
949 :
950 0 : ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3)
951 :
952 : !The dielectric stuff is performed in sequential mode.
953 : !Set mpi_enreg_diel accordingly
954 0 : call initmpi_seq(mpi_enreg_diel)
955 0 : call mpi_enreg_diel%distribfft%init_seq('c',ngfftdiel(2),ngfftdiel(3),'all')
956 0 : me_bandfft=xmpi_comm_rank(mpi_enreg%comm_bandfft)
957 :
958 0 : testocc=1
959 : !DEBUG
960 : !write(std_out,*)' susk : set testocc to 0 '
961 : !testocc=0
962 : !write(std_out,*)' susk : set extrap to 0 '
963 : !extrap=0
964 : !ENDDEBUG
965 :
966 : !Allocations, initializations
967 0 : ABI_MALLOC(rhoaug,(ndiel4,ndiel5,ndiel6))
968 0 : ABI_MALLOC(wfraug,(2,ndiel4,ndiel5,ndiel6))
969 0 : ABI_MALLOC(wfprod,(2,npwdiel))
970 0 : ABI_MALLOC(wfrspa,(2,ndiel4,ndiel5,ndiel6,nspinor,mband))
971 0 : ABI_MALLOC(dummy,(2,1))
972 0 : wfrspa(:,:,:,:,:,:)=zero
973 0 : ABI_MALLOC(treat_band,(nband_k))
974 0 : treat_band(:)=.true.
975 0 : isp1=isp;isp2=isp
976 0 : if (nspden_eff==2.and.nspinor==2) isp2=isp+1
977 :
978 : !BAND-FFT parallelism
979 0 : if (mpi_enreg%paral_kgb==1) then
980 0 : treat_band(:)=.false.
981 : ! We gather the wavefunctions treated by this proc in cg_local
982 0 : spaceComm=mpi_enreg%comm_band
983 0 : blocksize=mpi_enreg%nproc_band
984 0 : nbdblock=nband_k/blocksize
985 0 : ABI_MALLOC(sdispls,(blocksize))
986 0 : ABI_MALLOC(sdisplsloc,(blocksize))
987 0 : ABI_MALLOC(sendcounts,(blocksize))
988 0 : ABI_MALLOC(sendcountsloc,(blocksize))
989 0 : ABI_MALLOC(rdispls,(blocksize))
990 0 : ABI_MALLOC(rdisplsloc,(blocksize))
991 0 : ABI_MALLOC(recvcounts,(blocksize))
992 0 : ABI_MALLOC(recvcountsloc,(blocksize))
993 : ! First gather the kg_k in kg_k_gather_all
994 0 : npw_k_loc=npw_k_mpi
995 0 : call xmpi_allgather(npw_k_loc,recvcounts,spaceComm,ier)
996 0 : rdispls(1)=0
997 0 : do iproc=2,blocksize
998 0 : rdispls(iproc)=rdispls(iproc-1)+recvcounts(iproc-1)
999 : end do
1000 0 : ndatarecv=rdispls(blocksize)+recvcounts(blocksize)
1001 0 : ABI_MALLOC(kg_k_gather,(3,ndatarecv))
1002 0 : recvcountsloc(:)=recvcounts(:)*3
1003 0 : rdisplsloc(:)=rdispls(:)*3
1004 0 : call xmpi_allgatherv(kg_k_mpi,3*npw_k_loc,kg_k_gather,recvcountsloc(:),rdisplsloc,spaceComm,ier)
1005 0 : ABI_MALLOC(npw_per_proc,(mpi_enreg%nproc_fft))
1006 0 : ABI_MALLOC(rdispls_all,(mpi_enreg%nproc_fft))
1007 0 : spaceComm=mpi_enreg%comm_fft
1008 0 : call xmpi_allgather(ndatarecv,npw_per_proc,spaceComm,ier)
1009 0 : rdispls_all(1)=0
1010 0 : do iproc=2,mpi_enreg%nproc_fft
1011 0 : rdispls_all(iproc)=rdispls_all(iproc-1)+npw_per_proc(iproc-1)
1012 : end do
1013 0 : npw_tot=rdispls_all(mpi_enreg%nproc_fft)+npw_per_proc(mpi_enreg%nproc_fft)
1014 0 : ABI_MALLOC(kg_k_gather_all,(3,npw_tot))
1015 0 : call xmpi_allgatherv(kg_k_gather,3*ndatarecv,kg_k_gather_all,3*npw_per_proc(:),3*rdispls_all,spaceComm,ier)
1016 : ! At this point kg_k_gather_all contains all the kg
1017 : if(allocated(cwavef)) then
1018 : ABI_FREE(cwavef)
1019 : end if
1020 0 : ABI_MALLOC(cwavef,(2,npw_k_loc*nspinor*blocksize))
1021 0 : sizemax_per_proc=nband_k/(mpi_enreg%nproc_band*mpi_enreg%nproc_fft)+1
1022 0 : ABI_MALLOC(band_loc,(sizemax_per_proc))
1023 0 : ABI_MALLOC(cg_local,(2,sizemax_per_proc*npw_tot*nspinor))
1024 0 : iband_loc=0
1025 0 : do ibdblock=1,nbdblock
1026 : cwavef(:,1:npw_k_loc*nspinor*blocksize)=&
1027 0 : & cg_mpi(:,1+(ibdblock-1)*npw_k_loc*nspinor*blocksize+icg_mpi:ibdblock*npw_k_loc*nspinor*blocksize+icg_mpi)
1028 0 : sendcounts(:)=npw_k_loc
1029 0 : do iproc=1,blocksize
1030 0 : sdispls(iproc)=(iproc-1)*npw_k_loc
1031 : end do
1032 0 : ABI_MALLOC(cwavef_alltoall,(2,ndatarecv*nspinor))
1033 0 : recvcountsloc(:)=recvcounts(:)*2*nspinor
1034 0 : rdisplsloc(:)=rdispls(:)*2*nspinor
1035 0 : sendcountsloc(:)=sendcounts(:)*2*nspinor
1036 0 : sdisplsloc(:)=sdispls(:)*2*nspinor
1037 0 : call timab(547,1,tsec)
1038 0 : spaceComm=mpi_enreg%comm_band
1039 0 : call xmpi_alltoallv(cwavef,sendcountsloc,sdisplsloc,cwavef_alltoall,recvcountsloc,rdisplsloc,spaceComm,ier)
1040 0 : call timab(547,2,tsec)
1041 0 : ABI_MALLOC(cwavef_alltoall_gather,(2,npw_tot*nspinor))
1042 0 : blocksize=mpi_enreg%nproc_band
1043 0 : spaceComm=mpi_enreg%comm_fft
1044 : call xmpi_allgatherv(cwavef_alltoall,2*nspinor*ndatarecv,cwavef_alltoall_gather,&
1045 0 : & 2*nspinor*npw_per_proc,2*nspinor*rdispls_all,spaceComm,ier)
1046 0 : iproc_fft=modulo(ibdblock-1,mpi_enreg%nproc_fft)
1047 0 : if(mpi_enreg%me_fft==iproc_fft) then !All nproc_band procs of index me_fft will treat these bands
1048 0 : iband_loc=iband_loc+1
1049 0 : iband=1+mpi_enreg%me_band+mpi_enreg%nproc_band*mpi_enreg%me_fft+(iband_loc-1)*mpi_enreg%nproc_fft*mpi_enreg%nproc_band
1050 0 : treat_band(iband)=.true.
1051 0 : band_loc(iband_loc)=iband
1052 0 : cg_local(:,1+(iband_loc-1)*npw_tot*nspinor:iband_loc*npw_tot*nspinor)=cwavef_alltoall_gather(:,1:npw_tot*nspinor)
1053 : end if
1054 0 : ABI_FREE(cwavef_alltoall_gather)
1055 0 : ABI_FREE(cwavef_alltoall)
1056 : end do
1057 : ! On exit:
1058 : ! npw_tot will be npw
1059 : ! kg_k_gather_all will be kg_k
1060 : ! cg_local will be cg
1061 : ! icg will be zero
1062 0 : npw_k=>npw_tot
1063 0 : kg_k=>kg_k_gather_all(:,:)
1064 0 : cg=>cg_local(:,:)
1065 0 : icg=>icg_loc
1066 0 : call sphereboundary(gbound,istwf_k,kg_k,mgfftdiel,npw_k)
1067 0 : ABI_FREE(npw_per_proc)
1068 0 : ABI_FREE(rdispls_all)
1069 0 : ABI_FREE(sendcounts)
1070 0 : ABI_FREE(recvcounts)
1071 0 : ABI_FREE(sdispls)
1072 0 : ABI_FREE(rdispls)
1073 0 : ABI_FREE(sendcountsloc)
1074 0 : ABI_FREE(sdisplsloc)
1075 0 : ABI_FREE(recvcountsloc)
1076 0 : ABI_FREE(rdisplsloc)
1077 0 : ABI_FREE(kg_k_gather)
1078 0 : ABI_FREE(cwavef)
1079 : ! Because they will be summed over all procs, and arrive on input, rescale drhode and rhoextrap
1080 0 : if(occopt>=3)drhode(:,:,isp1:isp2)=drhode(:,:,isp1:isp2)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp)
1081 0 : if(extrap==1)rhoextrap(:,:,:,:)=rhoextrap(:,:,:,:)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp)
1082 0 : do i1=isp1,isp2
1083 0 : susmat(:,:,i1,:,i1)=susmat(:,:,i1,:,i1)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp)
1084 : end do
1085 :
1086 : ! No BAND-FFT parallelism
1087 : else ! use argument variables
1088 0 : cg=>cg_mpi
1089 0 : kg_k=>kg_k_mpi
1090 0 : npw_k=>npw_k_mpi
1091 0 : icg=>icg_mpi
1092 : end if
1093 0 : iband_loc=0
1094 :
1095 0 : call timab(751,2,tsec)
1096 0 : call timab(752,1,tsec)
1097 :
1098 : !Loop over bands to fft and store Fourier transform of wavefunction
1099 0 : ABI_MALLOC(cwavef,(2,npw_k))
1100 0 : do iband=1,nband_k
1101 0 : if(.not. treat_band(iband)) cycle ! I am not treating this band (only for the parallel case)
1102 0 : iband_loc=iband_loc+1
1103 :
1104 : ! Loop on spinorial components
1105 0 : do ispinor=1,nspinor
1106 0 : iwf=(ispinor-1)*npw_k+(iband_loc-1)*npw_k*nspinor+icg
1107 0 : jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1108 :
1109 : ! Obtain Fourier transform in fft box
1110 0 : tim_fourwf=8
1111 0 : cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf)
1112 : call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,&
1113 : & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,ndiel4,ndiel5,ndiel6,&
1114 0 : & 0,tim_fourwf,weight,weight)
1115 :
1116 0 : wfrspa(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:)
1117 :
1118 0 : if( (occopt>=3 .and. testocc==1) .or. extrap==1 )then
1119 : ! In the case of metallic occupation, or if the extrapolation
1120 : ! over higher bands is included, must compute the
1121 : ! Fourier transform of the density of each band, then
1122 : ! generate the part of the susceptibility matrix due
1123 : ! varying occupation numbers.
1124 :
1125 0 : weight=-two*occ_deavg(iband)*wtk(ikpt)/ucvol
1126 0 : do i3=1,ndiel3
1127 0 : do i2=1,ndiel2
1128 0 : do i1=1,ndiel1
1129 0 : wfraug(1,i1,i2,i3)=wfraug(1,i1,i2,i3)**2+wfraug(2,i1,i2,i3)**2
1130 0 : wfraug(2,i1,i2,i3)=zero
1131 : end do
1132 : end do
1133 : ! If extrapolation, accumulate density in real space
1134 0 : if(extrap==1.and.usepaw==0)then
1135 0 : do i2=1,ndiel2
1136 0 : do i1=1,ndiel1
1137 0 : rhoextrap(i1,i2,i3,ispinor)=rhoextrap(i1,i2,i3,ispinor)+weight*wfraug(1,i1,i2,i3)
1138 : end do
1139 : end do
1140 : end if
1141 : end do
1142 :
1143 : ! In case of PAW, add compensation charge contribution
1144 0 : if (usepaw==1.and.extrap==1.and.neglect_pawhat==0) then
1145 : call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,ispinor,ispinor,1,kg_diel,&
1146 : & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1147 : & ngfftdiel,npwdiel,nspinor,ntypat,1,&
1148 : & pawang,pawtab,ph3d_diel,typat,dummy,wfraug,&
1149 0 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1150 0 : rhoextrap(:,:,:,ispinor)=rhoextrap(:,:,:,ispinor)+weight*wfraug(1,:,:,:)
1151 : end if
1152 :
1153 : ! Performs the Fourier Transform of the density of the band,
1154 : ! and store it in wfprod
1155 0 : tim_fourwf=9
1156 : call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1157 : & 1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,&
1158 0 : & ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight)
1159 : ! In case of PAW, add compensation charge contribution if not already done
1160 0 : if (usepaw==1.and.extrap==0.and.neglect_pawhat==0) then
1161 : call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,&
1162 : & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1163 : & ngfftdiel,npwdiel,nspinor,ntypat,0,&
1164 : & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1165 0 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1166 : end if
1167 :
1168 : ! Perform now the summation of terms related to direct change of eigenvalues
1169 : ! or extrapolation over higher bands
1170 0 : wght1=zero ; wght2=zero
1171 0 : if(occopt>=3 .and. testocc==1) wght1=doccde(iband+bdtot_index)*wtk(ikpt)/ucvol
1172 0 : if(extrap==1) wght2=two*occ_deavg(iband)*wtk(ikpt)/ucvol
1173 0 : weight=wght1+wght2
1174 :
1175 0 : if (abs(weight)>tol12) then
1176 0 : do ipw2=1,npwdiel
1177 : ! Only fills lower half of the matrix (here, the susceptibility matrix)
1178 : ! Note that wfprod of the first index must behave like a density,
1179 : ! so that it is used as generated by fourwf, while wfprod of the
1180 : ! second index will be implicitely used to make a scalar product
1181 : ! with a potential change, meaning that its complex conjugate must be
1182 : ! used. This explains the following signs...
1183 0 : do ipw1=ipw2,npwdiel
1184 : susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1185 0 : & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1186 : susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1187 0 : & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1188 : end do
1189 : end do
1190 : end if
1191 :
1192 0 : if( occopt>=3 .and. testocc==1 .and. abs(wght1)>tol12) then
1193 : ! Accumulate product of band densities by their doccde, for the
1194 : ! computation of the effect of change of Fermi level.
1195 0 : do ipw=1,npwdiel
1196 0 : drhode(1,ipw,jsp)=drhode(1,ipw,jsp)+wfprod(1,ipw)*wght1
1197 0 : drhode(2,ipw,jsp)=drhode(2,ipw,jsp)+wfprod(2,ipw)*wght1
1198 : end do
1199 : ! Also accumulate weighted sum of doccde
1200 0 : sumdocc=sumdocc+wght1
1201 : end if
1202 :
1203 : ! End condition of metallic occupancies or extrapolation
1204 : end if
1205 :
1206 : ! End loop on spinorial components
1207 : end do
1208 : ! End loop on iband
1209 : end do
1210 :
1211 0 : call timab(752,2,tsec)
1212 0 : call timab(753,1,tsec)
1213 :
1214 0 : ABI_FREE(cwavef)
1215 :
1216 : !Stuff for parallelism (bands-FFT)
1217 0 : if(mpi_enreg%paral_kgb==1) then
1218 0 : call xmpi_sum(wfrspa,mpi_enreg%comm_bandfft,ier)
1219 0 : if(occopt>=3) then
1220 0 : call xmpi_sum(drhode(:,:,isp1:isp2),mpi_enreg%comm_bandfft,ier)
1221 : end if
1222 0 : if(extrap==1) then
1223 0 : call xmpi_sum(rhoextrap,mpi_enreg%comm_bandfft,ier)
1224 : end if
1225 0 : if(occopt>=3) then
1226 0 : call xmpi_sum(sumdocc,mpi_enreg%comm_bandfft,ier)
1227 : end if
1228 0 : ABI_MALLOC(susmat_mpi,(2,npwdiel,npwdiel))
1229 0 : do i1=isp1,isp2
1230 0 : susmat_mpi(:,:,:)=susmat(:,:,i1,:,i1)
1231 0 : call xmpi_sum(susmat_mpi,mpi_enreg%comm_bandfft,ier)
1232 0 : susmat(:,:,i1,:,i1)=susmat_mpi(:,:,:)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp)
1233 : end do
1234 0 : ABI_FREE(susmat_mpi)
1235 : end if
1236 0 : call timab(753,2,tsec)
1237 :
1238 : !-- Wavefunctions have been generated in real space ------------------------
1239 : !-- Now, compute product of wavefunctions for different bands --------------
1240 0 : call timab(754,1,tsec)
1241 : !if (occopt<3) then
1242 0 : tolocc=1.0d-3
1243 : !else
1244 : !tolocc=1.0d-8
1245 : !end if
1246 0 : iproc=-1
1247 :
1248 0 : if(nband_k>1)then
1249 0 : do ibd1=1,nband_k-1
1250 0 : do ibd2=ibd1+1,nband_k
1251 0 : iproc=iproc+1
1252 0 : if(modulo(iproc,mpi_enreg%nproc_fft*mpi_enreg%nproc_band) /= me_bandfft) cycle
1253 : ! If the occupation numbers are sufficiently different, or
1254 : ! if extrapolation is used and the corresponding factor is not zero,
1255 : ! then there is a contribution
1256 0 : occdiff=occ(ibd1+bdtot_index)-occ(ibd2+bdtot_index)
1257 0 : if( abs(occdiff)>tolocc .or. &
1258 : & ( extrap==1 .and. &
1259 : & ( abs(occ_deavg(ibd1)) + abs(occ_deavg(ibd2)) ) >tolocc ) &
1260 0 : & ) then
1261 :
1262 0 : eigdiff=eigen(ibd1+bdtot_index)-eigen(ibd2+bdtot_index)
1263 : ! DEBUG
1264 : ! write(std_out,*)' susk : contribution from bands',ibd1,ibd2
1265 : ! write(std_out,*)' occ diff =',occdiff
1266 : ! write(std_out,*)' eig diff =',eigdiff
1267 : ! ENDDEBUG
1268 :
1269 : ! Loop on spinorial components
1270 0 : do ispinor=1,nspinor
1271 0 : jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1272 :
1273 : ! Store the contribution in wfraug
1274 0 : do i3=1,ndiel3
1275 0 : do i2=1,ndiel2
1276 0 : do i1=1,ndiel1
1277 : wfraug(1,i1,i2,i3)=wfrspa(1,i1,i2,i3,ispinor,ibd1)*wfrspa(1,i1,i2,i3,ispinor,ibd2)&
1278 0 : & +wfrspa(2,i1,i2,i3,ispinor,ibd1)*wfrspa(2,i1,i2,i3,ispinor,ibd2)
1279 : wfraug(2,i1,i2,i3)=wfrspa(2,i1,i2,i3,ispinor,ibd1)*wfrspa(1,i1,i2,i3,ispinor,ibd2)&
1280 0 : & -wfrspa(1,i1,i2,i3,ispinor,ibd1)*wfrspa(2,i1,i2,i3,ispinor,ibd2)
1281 : end do
1282 : end do
1283 : end do
1284 :
1285 : ! Performs the Fourier Transform of the product, and store it in wfprod
1286 0 : tim_fourwf=19
1287 : call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1288 : & 1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,&
1289 0 : & ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight)
1290 :
1291 : ! In case of PAW, add compensation charge contribution
1292 0 : if (usepaw==1.and.neglect_pawhat==0) then
1293 : call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,&
1294 : & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1295 : & ngfftdiel,npwdiel,nspinor,ntypat,0,&
1296 : & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1297 0 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1298 : end if
1299 :
1300 : ! Perform now the summation
1301 0 : wght1=zero ; wght2=zero
1302 0 : if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol
1303 0 : if(extrap==1) wght2=(occ_deavg(ibd1)+occ_deavg(ibd2)) * two*wtk(ikpt)/ucvol
1304 0 : weight=wght1+wght2
1305 :
1306 : ! DEBUG
1307 : ! write(std_out,*)' weight =',weight
1308 : ! norm=zero
1309 : ! do ipw=1,npwdiel
1310 : ! norm=norm+wfprod(1,ipw)**2+wfprod(2,ipw)**2
1311 : ! end do
1312 : ! write(std_out,*)' norm in reciprocal space =',norm
1313 : ! ENDDEBUG
1314 :
1315 0 : if (abs(weight)>tol12) then
1316 0 : do ipw2=1,npwdiel
1317 : ! Only fills lower half of the matrix (here, the susceptibility matrix)
1318 : ! Note that wfprod of the first index must behave like a density,
1319 : ! so that it is used as generated by fourwf, while wfprod of the
1320 : ! second index will be implicitely used to make a scalar product
1321 : ! with a potential change, meaning that its complex conjugate must be
1322 : ! used. This explains the following signs...
1323 0 : do ipw1=ipw2,npwdiel
1324 : susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1325 0 : & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1326 : susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1327 0 : & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1328 : end do
1329 : end do
1330 : end if
1331 :
1332 : ! End loop on spinorial components
1333 : end do
1334 : ! End condition of different occupation numbers or extrapolation
1335 : end if
1336 : ! End internal loop over bands
1337 : end do
1338 : ! End external loop over bands
1339 : end do
1340 : ! End condition of having more than one band
1341 : end if
1342 :
1343 0 : call timab(754,2,tsec)
1344 0 : call timab(755,1,tsec)
1345 :
1346 0 : if(mpi_enreg%paral_kgb==1) then
1347 0 : ABI_MALLOC(susmat_mpi,(2,npwdiel,npwdiel))
1348 0 : do i1=isp1,isp2
1349 0 : susmat_mpi(:,:,:)=susmat(:,:,i1,:,i1)
1350 0 : call xmpi_sum(susmat_mpi,mpi_enreg%comm_bandfft,ier)
1351 0 : susmat(:,:,i1,:,i1)=susmat_mpi(:,:,:)
1352 : end do
1353 0 : ABI_FREE(susmat_mpi)
1354 0 : ABI_FREE(band_loc)
1355 0 : ABI_FREE(treat_band)
1356 0 : ABI_FREE(cg_local)
1357 0 : ABI_FREE(kg_k_gather_all)
1358 : end if
1359 :
1360 0 : call destroy_mpi_enreg(mpi_enreg_diel)
1361 0 : ABI_FREE(dummy)
1362 0 : ABI_FREE(rhoaug)
1363 0 : ABI_FREE(wfprod)
1364 0 : ABI_FREE(wfraug)
1365 0 : ABI_FREE(wfrspa)
1366 :
1367 0 : call timab(755,2,tsec)
1368 0 : call timab(750,2,tsec)
1369 :
1370 0 : end subroutine susk
1371 : !!***
1372 :
1373 : !!****f* m_suscep_stat/suskmm
1374 : !! NAME
1375 : !! suskmm
1376 : !!
1377 : !! FUNCTION
1378 : !! Compute the contribution of one k point to the susceptibility matrix
1379 : !! from input wavefunctions, band occupations, and k point wts.
1380 : !! Include the usual sum-over-state terms, but also the
1381 : !! corrections due to the change of the Fermi level in the metallic
1382 : !! case, as well as implicit sum over higher lying conduction
1383 : !! states, thanks to the closure relation (referred to as an extrapolation).
1384 : !!
1385 : !! This routine is similar to susk, but use blocking on wavefunctions
1386 : !! to decrease memory requirements, at the expense of CPU time.
1387 : !!
1388 : !! NOTES
1389 : !! There is still room for optimization !!
1390 : !!
1391 : !! INPUTS
1392 : !! atindx(natom)=index table for atoms
1393 : !! bdtot_index=index for the number of the band
1394 : !! cg(2,mcg)=wf in G space
1395 : !! cprj_k(natom,nspinor*nband_k)= wave functions projected with non-local projectors:
1396 : !! cprj_k=<p_i|Cnk> where p_i is a non-local projector.
1397 : !! doccde(mband*nkpt*nsppol)=derivative of occupancies wrt
1398 : !! the energy for each band and k point
1399 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
1400 : !! extrap: if==1, the closure relation (an extrapolation) must be used
1401 : !! gbound(2*mgfftdiel+8,2)=G sphere boundary for going from WF sphere to
1402 : !! medium size FFT grid
1403 : !! gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for going from medium size
1404 : !! FFT grid to small sphere.
1405 : !! gylmg_diel(npwdiel,lmax_diel**2,ntypat*usepaw)= -PAW only- Fourier transform of g_l(r).Y_ml(r) shape functions
1406 : !! for dielectric matrix
1407 : !! icg=index for cg
1408 : !! ikpt=number of the k point
1409 : !! isp=number of the current spin
1410 : !! istwf_k=input option parameter that describes the storage of wfs
1411 : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
1412 : !! kg_k(3,npw)=coordinates of planewaves in basis sphere.
1413 : !! lmax_diel=1+max. value of l angular momentum used for dielectric matrix
1414 : !! mband=maximum number of bands
1415 : !! mcg=dimension of cg
1416 : !! mgfftdiel=maximum size of 1D FFTs, for the computation of
1417 : !! the dielectric matrix
1418 : !! mpi_enreg=information about MPI parallelization
1419 : !! natom=number of atoms in cell
1420 : !! nband_k=number of bands at this k point for that spin polarization
1421 : !! ndiel4,ndiel5,ndiel6= FFT dimensions, modified to avoid cache trashing
1422 : !! neglect_pawhat=1 if PAW contribution from hat density (compensation charge)
1423 : !! has to be neglected (to be used when only an estimation of
1424 : !! suscep. matrix has to be evaluated, i.e. for SCF precondictioning)
1425 : !! ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix,
1426 : !! see ~abinit/doc/variables/vargs.htm#ngfft
1427 : !! nkpt=number of k points
1428 : !! npwdiel=third and fifth dimension of the susmat array.
1429 : !! npw_k=number of plane waves at this k point
1430 : !! nspden=number of spin-density components
1431 : !! nspden_eff=number of spin-density components actually computed in sussceptibility
1432 : !! nspinor=number of spinorial components of the wavefunctions
1433 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1434 : !! ntypat=number of types of atoms in unit cell.
1435 : !! occ(mband*nkpt*nsppol)=
1436 : !! occupation numbers for each band (usually 2.0) at each k point
1437 : !! occopt=option for occupancies
1438 : !! occ_deavg(mband)=factor for extrapolation (occup. divided by an energy gap)
1439 : !! pawang <type(pawang_type)>=paw angular mesh and related data
1440 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
1441 : !! ph3d_diel(2,npwdiel,natom*usepaw)=3-dim structure factors, for each atom and plane wave, for dielectric matrix
1442 : !! typat(natom)=type (integer) for each atom
1443 : !! ucvol=unit cell volume (Bohr**3)
1444 : !! usepaw=flag for PAW
1445 : !! wtk(nkpt)=k point weights (they sum to 1.0)
1446 : !!
1447 : !! OUTPUT
1448 : !! (see side effects)
1449 : !!
1450 : !! SIDE EFFECTS
1451 : !! drhode(2,npwdiel,nspden_eff)=weighted density, needed to compute the
1452 : !! effect of change of fermi energy
1453 : !! rhoextrap(ndiel4,ndiel5,ndiel6,nspinor)=density-like array, needed for the
1454 : !! extrapolation procedure.
1455 : !! sumdocc=sum of weighted occupation numbers, needed to compute the
1456 : !! effect of change of fermi energy
1457 : !! susmat(2,npwdiel,nspden,npwdiel,nspden)=
1458 : !! the susceptibility (or density-density response) matrix in reciprocal space
1459 : !!
1460 : !! SOURCE
1461 :
1462 322 : subroutine suskmm(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,&
1463 322 : & gbound_diel,gylmg_diel,icg,ikpt,isp,istwf_k,kg_diel,kg_k,&
1464 : & lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,&
1465 : & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,&
1466 322 : & npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,&
1467 322 : & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,&
1468 322 : & susmat,typat,ucvol,usepaw,wtk)
1469 :
1470 : !Arguments ------------------------------------
1471 : !scalars
1472 : integer,intent(in) :: bdtot_index,extrap,icg,ikpt,isp,istwf_k,lmax_diel,mband,mcg
1473 : integer,intent(in) :: mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat
1474 : integer,intent(in) :: nkpt,npw_k,npwdiel,nspden,nspden_eff,nspinor
1475 : integer,intent(in) :: nsppol,ntypat,occopt,usepaw
1476 : real(dp),intent(in) :: ucvol
1477 : real(dp),intent(inout) :: sumdocc
1478 : type(MPI_type),intent(in) :: mpi_enreg
1479 : type(pawang_type),intent(in) :: pawang
1480 : !arrays
1481 : integer,intent(in) :: atindx(natom),gbound(2*mgfftdiel+8,2)
1482 : integer,intent(in) :: gbound_diel(2*mgfftdiel+8,2)
1483 : integer,intent(in) :: kg_diel(3,npwdiel),kg_k(3,npw_k),ngfftdiel(18)
1484 : integer,intent(in) :: typat(natom)
1485 : real(dp),intent(in) :: cg(2,mcg),doccde(mband*nkpt*nsppol)
1486 : real(dp),intent(in) :: eigen(mband*nkpt*nsppol)
1487 : real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat*usepaw)
1488 : real(dp),intent(in) :: occ(mband*nkpt*nsppol),occ_deavg(mband)
1489 : real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom*usepaw),wtk(nkpt)
1490 : real(dp),intent(inout) :: drhode(2,npwdiel,nspden_eff)
1491 : real(dp),intent(inout) :: rhoextrap(ndiel4,ndiel5,ndiel6,nspinor)
1492 : real(dp),intent(inout) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
1493 : type(pawcprj_type) :: cprj_k(natom,nspinor*nband_k*usepaw)
1494 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
1495 :
1496 : !Local variables-------------------------------
1497 : !scalars
1498 : integer :: comm_fft,i1,i2,i3,iband,iband_shift,iband_shift2,ibd1,ibd2,ibdshft1,ibdshft2
1499 : integer :: iblk1,iblk2,ipw,ipw1,ipw2,ispinor,iwf,jsp,mblk
1500 : integer :: nblk,nbnd_current,nbnd_in_blk,nbnd_in_blk1,ndiel1,ndiel2,ndiel3
1501 : integer :: testocc,tim_fourwf
1502 : real(dp) :: eigdiff,occdiff,tolocc,weight,wght1,wght2
1503 : character(len=500) :: message
1504 322 : type(MPI_type) :: mpi_enreg_diel
1505 : !arrays
1506 : real(dp) :: tsec(2)
1507 322 : real(dp),allocatable :: cwavef(:,:),dummy(:,:),rhoaug(:,:,:),wfprod(:,:)
1508 322 : real(dp),allocatable :: wfraug(:,:,:,:),wfrspa1(:,:,:,:,:,:)
1509 322 : real(dp),allocatable :: wfrspa2(:,:,:,:,:,:)
1510 :
1511 : ! *************************************************************************
1512 :
1513 322 : call timab(760,1,tsec)
1514 322 : call timab(761,1,tsec)
1515 :
1516 : !Allocations, initializations
1517 322 : ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3)
1518 322 : testocc=1
1519 1610 : ABI_MALLOC(rhoaug,(ndiel4,ndiel5,ndiel6))
1520 1610 : ABI_MALLOC(wfraug,(2,ndiel4,ndiel5,ndiel6))
1521 966 : ABI_MALLOC(wfprod,(2,npwdiel))
1522 322 : ABI_MALLOC(dummy,(2,1))
1523 :
1524 : !The dielectric stuff is performed in sequential mode.
1525 : !Set mpi_enreg_diel accordingly
1526 322 : call initmpi_seq(mpi_enreg_diel)
1527 322 : call mpi_enreg_diel%distribfft%init_seq('c',ngfftdiel(2),ngfftdiel(3),'all')
1528 :
1529 322 : comm_fft=mpi_enreg%comm_fft
1530 :
1531 : !Prepare the blocking : compute the number of blocks,
1532 : !the number of bands in each normal block,
1533 : !and the number in the first one, usually smaller.
1534 :
1535 : !Consider that if the number of bands is large, there are at most 8 blocks
1536 322 : nbnd_in_blk=0
1537 322 : if(nband_k>=48)then
1538 0 : mblk=8
1539 0 : nbnd_in_blk=(nband_k-1)/mblk+1
1540 : ! If the number of bands is medium, place 6 bands per block
1541 322 : else if(nband_k>=12)then
1542 : nbnd_in_blk=6
1543 : ! Otherwise, must have at least 2 blocks
1544 174 : else if(nband_k>=2)then
1545 174 : mblk=2
1546 174 : nbnd_in_blk=(nband_k-1)/mblk+1
1547 : else
1548 : write(message, '(a,a,a,i2,a,a,a)')&
1549 0 : & ' The number of bands must be larger or equal to 2, in suskmm.',ch10,&
1550 0 : & ' It is equal to ',nband_k,'.',ch10,&
1551 0 : & ' Action : choose another preconditioner.'
1552 0 : ABI_ERROR(message)
1553 : end if
1554 :
1555 : !Compute the effective number of blocks, and the number of bands in
1556 : !the first block.
1557 322 : nblk=(nband_k-1)/nbnd_in_blk+1
1558 322 : nbnd_in_blk1=nband_k-(nblk-1)*nbnd_in_blk
1559 :
1560 : !DEBUG
1561 : !write(std_out,*)' suskmm : nband_k,nblk,nbnd_in_blk,nbnd_in_blk1 '
1562 : !write(std_out,*)nband_k,nblk,nbnd_in_blk,nbnd_in_blk1
1563 : !stop
1564 : !ENDDEBUG
1565 :
1566 : !wfrspa1 will contain the wavefunctions of the slow sampling (iblk1)
1567 2254 : ABI_MALLOC(wfrspa1,(2,ndiel4,ndiel5,ndiel6,nspinor,nbnd_in_blk))
1568 : !wfrspa2 will contain the wavefunctions of the rapid sampling (iblk2)
1569 1932 : ABI_MALLOC(wfrspa2,(2,ndiel4,ndiel5,ndiel6,nspinor,nbnd_in_blk))
1570 :
1571 966 : ABI_MALLOC(cwavef,(2,npw_k))
1572 :
1573 322 : call timab(761,2,tsec)
1574 :
1575 : !First loop over blocks
1576 1146 : do iblk1=1,nblk
1577 :
1578 824 : call timab(762,1,tsec)
1579 :
1580 : ! Initialisation
1581 824 : if(iblk1==1)then
1582 :
1583 322 : nbnd_current=nbnd_in_blk1
1584 322 : iband_shift=0
1585 : ! Loop over bands to fft and store Fourier transform of wavefunction
1586 1518 : do iband=1,nbnd_current
1587 : ! Loop on spinorial components
1588 2714 : do ispinor=1,nspinor
1589 1196 : iwf=(ispinor-1)*npw_k+(iband-1)*npw_k*nspinor+icg
1590 : ! Obtain Fourier transform in fft box
1591 1196 : tim_fourwf=21
1592 1650728 : cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf)
1593 : call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,&
1594 : & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,ndiel4,ndiel5,ndiel6,&
1595 1196 : & 0,tim_fourwf,weight,weight)
1596 12394496 : wfrspa1(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:)
1597 : end do
1598 : end do
1599 :
1600 : else
1601 :
1602 : ! The Fourier transform of wavefunctions have already been obtained
1603 502 : nbnd_current=nbnd_in_blk
1604 502 : iband_shift=nbnd_in_blk1+(iblk1-2)*nbnd_in_blk
1605 :
1606 : end if
1607 :
1608 : ! Loop over bands of this block, to generate band-diagonal
1609 4800 : do iband=1,nbnd_current
1610 :
1611 : ! Loop on spinorial components
1612 8776 : do ispinor=1,nspinor
1613 3976 : jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1614 :
1615 7952 : if( (occopt>=3 .and. testocc==1) .or. extrap==1 )then
1616 : ! In the case of metallic occupation, or if the extrapolation
1617 : ! over higher bands is included, must compute the
1618 : ! Fourier transform of the density of each band, then
1619 : ! generate the part of the susceptibility matrix due
1620 : ! varying occupation numbers.
1621 3956 : weight=-two*occ_deavg(iband+iband_shift)*wtk(ikpt)/ucvol
1622 148632 : do i3=1,ndiel3
1623 1371076 : do i2=1,ndiel2
1624 11902916 : do i1=1,ndiel1
1625 : wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,iband)**2&
1626 10531840 : & +wfrspa1(2,i1,i2,i3,ispinor,iband)**2
1627 11758240 : wfraug(2,i1,i2,i3)=zero
1628 : end do
1629 : end do
1630 : ! If extrapolation, accumulate density in real space
1631 148632 : if(extrap==1.and.usepaw==0)then
1632 1371076 : do i2=1,ndiel2
1633 11902916 : do i1=1,ndiel1
1634 11758240 : rhoextrap(i1,i2,i3,ispinor)=rhoextrap(i1,i2,i3,ispinor)+weight*wfraug(1,i1,i2,i3)
1635 : end do
1636 : end do
1637 : end if
1638 : end do
1639 :
1640 : ! In case of PAW, add compensation charge contribution
1641 3956 : if (usepaw==1.and.extrap==1.and.neglect_pawhat==0) then
1642 : call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,ispinor,ispinor,1,kg_diel,&
1643 : & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1644 : & ngfftdiel,npwdiel,nspinor,ntypat,1,&
1645 : & pawang,pawtab,ph3d_diel,typat,dummy,wfraug,&
1646 0 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1647 0 : rhoextrap(:,:,:,ispinor)=rhoextrap(:,:,:,ispinor)+weight*wfraug(1,:,:,:)
1648 : end if
1649 :
1650 : ! Performs the Fourier Transform of the density of the band,
1651 : ! and store it in wfprod
1652 3956 : tim_fourwf=31
1653 : call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1654 : & 1,kg_diel,kg_diel,&
1655 3956 : & mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight)
1656 : ! In case of PAW, add compensation charge contribution if not already done
1657 3956 : if (usepaw==1.and.extrap==0.and.neglect_pawhat==0) then
1658 : call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,1,1,1,kg_diel,&
1659 : & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1660 : & ngfftdiel,npwdiel,nspinor,ntypat,0,&
1661 : & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1662 0 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1663 : end if
1664 :
1665 : ! Perform now the summation of terms related to direct change of eigenvalues
1666 : ! or extrapolation over higher bands
1667 3956 : wght1=zero ; wght2=zero
1668 3956 : if(occopt>=3 .and. testocc==1) wght1=doccde(iband+iband_shift+bdtot_index)*wtk(ikpt)/ucvol
1669 3956 : if(extrap==1) wght2=two*occ_deavg(iband+iband_shift)*wtk(ikpt)/ucvol
1670 3956 : weight=wght1+wght2
1671 :
1672 3956 : if (abs(weight)>tol12) then
1673 534560 : do ipw2=1,npwdiel
1674 : ! Only fills lower half of the matrix (here, the susceptibility matrix)
1675 : ! Note that wfprod of the first index must behave like a density,
1676 : ! so that it is used as generated by fourwf, while wfprod of the
1677 : ! second index will be implicitely used to make a scalar product
1678 : ! with a potential change, meaning that its complex conjugate must be
1679 : ! used. This explains the following signs...
1680 58617240 : do ipw1=ipw2,npwdiel
1681 : susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1682 58082680 : & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1683 : susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1684 58613559 : & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1685 : end do
1686 : end do
1687 : end if
1688 :
1689 3956 : if( occopt>=3 .and. testocc==1 .and. abs(wght1)>tol12) then
1690 : ! Accumulate product of band densities by their doccde, for the
1691 : ! computation of the effect of change of Fermi level.
1692 222102 : do ipw=1,npwdiel
1693 219926 : drhode(1,ipw,jsp)=drhode(1,ipw,jsp)+wfprod(1,ipw)*wght1
1694 222102 : drhode(2,ipw,jsp)=drhode(2,ipw,jsp)+wfprod(2,ipw)*wght1
1695 : end do
1696 : ! Also accumulate weighted sum of doccde
1697 2176 : sumdocc=sumdocc+wght1
1698 : end if
1699 :
1700 : ! End condition of metallic occupancies or extrapolation
1701 : end if
1702 :
1703 : ! End loop on spinorial components
1704 : end do
1705 : ! End loop on iband
1706 : end do
1707 :
1708 824 : call timab(762,2,tsec)
1709 :
1710 : ! -- Compute now off-band-diagonal terms ------------------------------------
1711 : ! -- Compute product of wavefunctions for different bands, inside the block -
1712 :
1713 824 : call timab(763,1,tsec)
1714 :
1715 : ! if (occopt<3) then
1716 824 : tolocc=1.0d-3
1717 : ! else
1718 : ! tolocc=1.0d-8
1719 : ! end if
1720 :
1721 824 : if(nbnd_current>1)then
1722 3976 : do ibd1=1,nbnd_current-1
1723 3152 : ibdshft1=ibd1+iband_shift
1724 12304 : do ibd2=ibd1+1,nbnd_current
1725 8328 : ibdshft2=ibd2+iband_shift
1726 :
1727 : ! If the occupation numbers are sufficiently different, or
1728 : ! if extrapolation is used and the corresponding factor is not zero,
1729 : ! then there is a contribution
1730 8328 : occdiff=occ(ibdshft1+bdtot_index)-occ(ibdshft2+bdtot_index)
1731 8328 : if( abs(occdiff)>tolocc .or. &
1732 : & ( extrap==1 .and. &
1733 : & ( abs(occ_deavg(ibdshft1)) + abs(occ_deavg(ibdshft2)) ) >tolocc ) &
1734 3152 : & ) then
1735 :
1736 7589 : eigdiff=eigen(ibdshft1+bdtot_index) - eigen(ibdshft2+bdtot_index)
1737 :
1738 : ! Loop on spinorial components
1739 15178 : do ispinor=1,nspinor
1740 7589 : jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1741 :
1742 : ! Store the contribution in wfraug
1743 291394 : do i3=1,ndiel3
1744 2711958 : do i2=1,ndiel2
1745 23602921 : do i1=1,ndiel1
1746 : wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa1(1,i1,i2,i3,ispinor,ibd2)&
1747 20898552 : & +wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa1(2,i1,i2,i3,ispinor,ibd2)
1748 : wfraug(2,i1,i2,i3)=wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa1(1,i1,i2,i3,ispinor,ibd2)&
1749 23319116 : & -wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa1(2,i1,i2,i3,ispinor,ibd2)
1750 : end do
1751 : end do
1752 : end do
1753 :
1754 : ! Performs the Fourier Transform of the product, and store it in wfprod
1755 7589 : tim_fourwf=32
1756 : call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1757 : & 1,kg_diel,kg_diel, mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,&
1758 7589 : & ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight)
1759 :
1760 : ! In case of PAW, add compensation charge contribution
1761 7589 : if (usepaw==1.and.neglect_pawhat==0) then
1762 : call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,&
1763 : & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1764 : & ngfftdiel,npwdiel,nspinor,ntypat,0,&
1765 : & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1766 0 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1767 : end if
1768 :
1769 : ! Perform now the summation
1770 7589 : wght1=zero ; wght2=zero
1771 7589 : if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol
1772 7589 : if(extrap==1) wght2=(occ_deavg(ibdshft1)+occ_deavg(ibdshft2)) * two*wtk(ikpt)/ucvol
1773 7589 : weight=wght1+wght2
1774 :
1775 15917 : if (abs(weight)>tol12) then
1776 1136872 : do ipw2=1,npwdiel
1777 : ! Only fills lower half of the matrix (here, the susceptibility matrix)
1778 : ! Note that wfprod of the first index must behave like a density,
1779 : ! so that it is used as generated by fourwf, while wfprod of the
1780 : ! second index will be implicitely used to make a scalar product
1781 : ! with a potential change, meaning that its complex conjugate must be
1782 : ! used. This explains the following signs...
1783 128619848 : do ipw1=ipw2,npwdiel
1784 : susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1785 127482976 : & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1786 : susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1787 128612259 : & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1788 : end do
1789 : end do
1790 : end if
1791 :
1792 : ! End loop on spinorial components
1793 : end do
1794 : ! End condition of different occupation numbers or extrapolation
1795 : end if
1796 : ! End internal loop over bands
1797 : end do
1798 : ! End external loop over bands
1799 : end do
1800 : ! End condition of having more than one band
1801 : end if
1802 :
1803 : ! Loop on secondary block, with fast varying index, in decreasing order.
1804 824 : if(iblk1/=nblk)then
1805 1228 : do iblk2=nblk,iblk1+1,-1
1806 726 : iband_shift2=nbnd_in_blk1+(iblk2-2)*nbnd_in_blk
1807 :
1808 : ! Loop over bands to fft and store Fourier transform of wavefunction
1809 726 : iband_shift2=nbnd_in_blk1+(iblk2-2)*nbnd_in_blk
1810 4850 : do iband=1,nbnd_in_blk
1811 : ! Loop on spinorial components
1812 8974 : do ispinor=1,nspinor
1813 4124 : iwf=(ispinor-1)*npw_k+(iband+iband_shift2-1)*npw_k*nspinor+icg
1814 :
1815 : ! Obtain Fourier transform in fft box
1816 4124 : tim_fourwf=22
1817 4975922 : cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf)
1818 : call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,&
1819 : & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,&
1820 4124 : & ndiel4,ndiel5,ndiel6,0,tim_fourwf,weight,weight)
1821 41682776 : wfrspa2(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:)
1822 : end do
1823 : end do
1824 :
1825 4270 : do ibd1=1,nbnd_current
1826 3042 : ibdshft1=ibd1+iband_shift
1827 21228 : do ibd2=1,nbnd_in_blk
1828 17460 : ibdshft2=ibd2+iband_shift2
1829 :
1830 : ! If the occupation numbers are sufficiently different, or
1831 : ! if extrapolation is used and the corresponding factor is not zero,
1832 : ! then there is a contribution
1833 17460 : occdiff=occ(ibdshft1+bdtot_index)-occ(ibdshft2+bdtot_index)
1834 17460 : if( abs(occdiff)>tolocc .or. &
1835 : & ( extrap==1 .and. &
1836 : & ( abs(occ_deavg(ibdshft1)) + abs(occ_deavg(ibdshft2)) ) >tolocc ) &
1837 3042 : & ) then
1838 :
1839 16481 : eigdiff=eigen(ibdshft1+bdtot_index) - eigen(ibdshft2+bdtot_index)
1840 :
1841 : ! Loop on spinorial components
1842 32962 : do ispinor=1,nspinor
1843 16481 : jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1844 :
1845 : ! Store the contribution in wfraug
1846 608213 : do i3=1,ndiel3
1847 5696749 : do i2=1,ndiel2
1848 50025916 : do i1=1,ndiel1
1849 : wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa2(1,i1,i2,i3,ispinor,ibd2)&
1850 44345648 : & +wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa2(2,i1,i2,i3,ispinor,ibd2)
1851 : wfraug(2,i1,i2,i3)=wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa2(1,i1,i2,i3,ispinor,ibd2)&
1852 49434184 : & -wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa2(2,i1,i2,i3,ispinor,ibd2)
1853 : end do
1854 : end do
1855 : end do
1856 :
1857 : ! Performs the Fourier Transform of the product, and store it in wfprod
1858 16481 : tim_fourwf=32
1859 : call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1860 : & 1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,&
1861 16481 : & ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight)
1862 :
1863 : ! In case of PAW, add compensation charge contribution
1864 16481 : if (usepaw==1.and.neglect_pawhat==0) then
1865 : call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibdshft2,ispinor,ispinor,1,kg_diel,&
1866 : & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1867 : & ngfftdiel,npwdiel,nspinor,ntypat,0,&
1868 : & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1869 0 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1870 : end if
1871 :
1872 : ! Perform now the summation
1873 16481 : wght1=zero ; wght2=zero
1874 16481 : if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol
1875 16481 : if(extrap==1) wght2=(occ_deavg(ibdshft1)+occ_deavg(ibdshft2)) * two*wtk(ikpt)/ucvol
1876 16481 : weight=wght1+wght2
1877 :
1878 33941 : if (abs(weight)>tol12) then
1879 2417690 : do ipw2=1,npwdiel
1880 : ! Only fills lower half of the matrix (here, the susceptibility matrix)
1881 : ! Note that wfprod of the first index must behave like a density,
1882 : ! so that it is used as generated by fourwf, while wfprod of the
1883 : ! second index will be implicitely used to make a scalar product
1884 : ! with a potential change, meaning that its complex conjugate must be
1885 : ! used. This explains the following signs...
1886 284800175 : do ipw1=ipw2,npwdiel
1887 : susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1888 282382485 : & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1889 : susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1890 284783694 : & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1891 : end do
1892 : end do
1893 : end if
1894 :
1895 : ! End loop on spinorial components
1896 : end do
1897 : ! End condition of different occupation numbers or extrapolation
1898 : end if
1899 : ! End internal loop over bands
1900 : end do
1901 : ! End external loop over bands
1902 : end do
1903 : ! End loop on bloks
1904 : end do
1905 :
1906 : ! Finish the loop on blok with iblk2=iblk1+1, so can use the
1907 : ! FFTd wavefunctions for the next iblk1.
1908 3282 : do iband=1,nbnd_in_blk
1909 28652370 : wfrspa1(:,:,:,:,1:nspinor,iband)=wfrspa2(:,:,:,:,1:nspinor,iband)
1910 : end do
1911 :
1912 : ! End condition of iblk1/=nblk
1913 : end if
1914 :
1915 1146 : call timab(763,2,tsec)
1916 :
1917 : ! End loop on iblk1
1918 : end do
1919 :
1920 : !DEBUG
1921 : !write(std_out,*)' suskmm : exit '
1922 : !do ipw1=1,npwdiel
1923 : !write(std_out,*)ipw1,susmat(1,ipw1,1,ipw1,1),susmat(2,ipw1,1,ipw1,1)
1924 : !end do
1925 : !write(std_out,*)' suskmm : end of susmat '
1926 : !stop
1927 : !ENDDEBUG
1928 :
1929 322 : call destroy_mpi_enreg(mpi_enreg_diel)
1930 322 : ABI_FREE(cwavef)
1931 322 : ABI_FREE(dummy)
1932 322 : ABI_FREE(rhoaug)
1933 322 : ABI_FREE(wfprod)
1934 322 : ABI_FREE(wfraug)
1935 322 : ABI_FREE(wfrspa1)
1936 322 : ABI_FREE(wfrspa2)
1937 :
1938 322 : call timab(760,2,tsec)
1939 :
1940 322 : end subroutine suskmm
1941 : !!***
1942 :
1943 : end module m_suscep_stat
1944 : !!***
|
