Branch data 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 : : !!***
|
