Line data Source code
1 : !!****m* ABINIT/m_paw_overlap
2 : !! NAME
3 : !! m_paw_overlap
4 : !!
5 : !! FUNCTION
6 : !! This module contains several routines used to compute the overlap between
7 : !! 2 wave-functions (PAW only), and associated tools.
8 : !! Mainly used in Berry phase formalism.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2018-2026 ABINIT group (JWZ,TRangel,BA,FJ,PHermet)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : MODULE m_paw_overlap
25 :
26 : use defs_basis
27 : use m_errors
28 : use m_abicore
29 : use m_xmpi
30 :
31 : use defs_abitypes, only : MPI_type
32 : use m_special_funcs, only : sbf8
33 : use m_efield, only : efield_type
34 : use m_pawang, only : pawang_type
35 : use m_pawcprj, only : pawcprj_type
36 : use m_pawrad, only : pawrad_type,simp_gen
37 : use m_pawtab, only : pawtab_type
38 : use m_paw_sphharm, only : initylmr
39 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free, pawcprj_getdim
40 :
41 : implicit none
42 :
43 : private
44 :
45 : !public procedures.
46 : public :: overlap_k1k2_paw
47 : public :: smatrix_pawinit
48 : public :: smatrix_k_paw
49 : public :: qijb_kk
50 : public :: expibi
51 :
52 : CONTAINS !========================================================================================
53 : !!***
54 :
55 : !----------------------------------------------------------------------
56 :
57 : !!****f* m_paw_overlap/overlap_k1k2_paw
58 : !! NAME
59 : !! overlap_k1k2_paw
60 : !!
61 : !! FUNCTION
62 : !! compute PAW overlap between two k points,
63 : !! similar to smatrix_k_paw.F90 but more generic
64 : !!
65 : !! INPUTS
66 : !! cprj_k1 (pawcprj_type) :: cprj for occupied bands at point k1
67 : !! cprj_k2 :: cprj for occupied bands at point k2
68 : !! dk(3) :: vector k2 - k1
69 : !! gprimd(3,3)=dimensioned primitive translations of reciprocal lattice
70 : !! lmn2max :: lmnmax*(lmnmax+1)/2
71 : !! lmnsize(ntypat) :: lmnsize for each atom type
72 : !! mband :: number of bands
73 : !! natom=number of atoms in unit cell
74 : !! nspinor :: number of spinors (1 or 2)
75 : !! ntypat=number of types of atoms in unit cell
76 : !! pawang <type(pawang_type)>=paw angular mesh and related data
77 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
78 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
79 : !! typat=typat(natom) list of atom types
80 : !! xred(natom,3) :: locations of atoms in cell
81 : !!
82 : !! OUTPUT
83 : !! k1k2_paw(2,mband,mband) :: array of the on-site PAW parts of the overlaps between Bloch states at points
84 : !! k1 and k2, for the various pairs of bands, that is, the on-site part of
85 : !! <u_nk1|u_mk2>
86 : !!
87 : !! SIDE EFFECTS
88 : !!
89 : !! NOTES
90 : !! This routine assumes that the cprj are not explicitly ordered by
91 : !! atom type.
92 : !!
93 : !! SOURCE
94 :
95 0 : subroutine overlap_k1k2_paw(cprj_k1,cprj_k2,dk,gprimd,k1k2_paw,lmn2max,lmnsize,&
96 0 : & natom,nband,nband_occ,nspinor,ntypat,pawang,pawrad,pawtab,typat,xred)
97 :
98 : !Arguments---------------------------
99 : !scalars
100 : integer,intent(in) :: lmn2max,natom,nband,nband_occ,nspinor,ntypat
101 : type(pawang_type),intent(in) :: pawang
102 : type(pawcprj_type),intent(in) :: cprj_k1(natom,nband),cprj_k2(natom,nband)
103 :
104 : !arrays
105 : integer,intent(in) :: lmnsize(ntypat),typat(natom)
106 : real(dp),intent(in) :: dk(3),gprimd(3,3),xred(natom,3)
107 : real(dp),intent(out) :: k1k2_paw(2,nband_occ,nband_occ)
108 : type(pawrad_type),intent(in) :: pawrad(ntypat)
109 : type(pawtab_type),intent(in) :: pawtab(ntypat)
110 :
111 : !Local variables---------------------------
112 : !scalars
113 : integer :: iatom,iband,ibs,ilmn,ispinor,itypat
114 : integer :: jband,jbs,jlmn,klmn
115 : complex(dp) :: cpk1,cpk2,cterm,paw_onsite
116 :
117 : ! arrays
118 0 : real(dp),allocatable :: calc_expibi(:,:),calc_qijb(:,:,:)
119 : ! *************************************************************************
120 :
121 : !initialize k1k2_paw output variable
122 0 : k1k2_paw(:,:,:) = zero
123 :
124 : ! obtain the atomic phase factors for the input k vector shift
125 0 : ABI_MALLOC(calc_expibi,(2,natom))
126 0 : call expibi(calc_expibi,dk,natom,xred)
127 :
128 : ! obtain the onsite PAW terms for the input k vector shift
129 0 : ABI_MALLOC(calc_qijb,(2,lmn2max,natom))
130 0 : call qijb_kk(calc_qijb,dk,calc_expibi,gprimd,lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat)
131 0 : ABI_FREE(calc_expibi)
132 :
133 0 : do iatom = 1, natom
134 0 : itypat = typat(iatom)
135 :
136 0 : do ilmn=1,lmnsize(itypat)
137 0 : do jlmn=1,lmnsize(itypat)
138 0 : klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn)
139 0 : paw_onsite = cmplx(calc_qijb(1,klmn,iatom),calc_qijb(2,klmn,iatom))
140 0 : do ispinor = 1, nspinor
141 0 : do iband = 1, nband_occ
142 0 : do jband = 1, nband_occ
143 0 : ibs = nspinor*(iband-1) + ispinor
144 0 : jbs = nspinor*(jband-1) + ispinor
145 0 : cpk1=cmplx(cprj_k1(iatom,ibs)%cp(1,ilmn),cprj_k1(iatom,ibs)%cp(2,ilmn))
146 0 : cpk2=cmplx(cprj_k2(iatom,jbs)%cp(1,jlmn),cprj_k2(iatom,jbs)%cp(2,jlmn))
147 0 : cterm = conjg(cpk1)*paw_onsite*cpk2
148 0 : k1k2_paw(1,iband,jband) = k1k2_paw(1,iband,jband)+real(cterm)
149 0 : k1k2_paw(2,iband,jband) = k1k2_paw(2,iband,jband)+aimag(cterm)
150 : end do ! end loop over jband
151 : end do ! end loop over iband
152 : end do ! end loop over ispinor
153 : end do ! end loop over ilmn
154 : end do ! end loop over jlmn
155 :
156 : end do ! end loop over atoms
157 :
158 0 : ABI_FREE(calc_qijb)
159 :
160 0 : end subroutine overlap_k1k2_paw
161 : !!***
162 :
163 : !----------------------------------------------------------------------
164 :
165 : !!****f* m_paw_overlap/smatrix_pawinit
166 : !! NAME
167 : !! smatrix_pawinit
168 : !!
169 : !! FUNCTION
170 : !! Routine which computes paw part of the overlap used to compute LMWF wannier
171 : !! functions and berryphase
172 : !!
173 : !! INPUTS
174 : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
175 : !! cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
176 : !! and each |p_lmn> non-local projector
177 : !! dimcprj(natom)=array of dimensions of array cprj (not ordered)
178 : !! g1(3)= reciprocal vector to put k1+b inside the BZ. bb=k2-k1=b-G1
179 : !! ("b" is the true b, so we have to correct bb with G1).
180 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
181 : !! ikpt1(3)=cartesian coordinates of k1
182 : !! ikpt2(3)=cartesian coordinates of k2
183 : !! isppol = spin polarization
184 : !! mband=maximum number of bands
185 : !! mkmem =number of k points treated by this node.
186 : !! mpi_enreg=information about MPI parallelization
187 : !! natom=number of atoms in cell.
188 : !! nkpt=number of k points.
189 : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
190 : !! nsppol=1 for unpolarized, 2 for spin-polarized
191 : !! ntypat=number of types of atoms in unit cell.
192 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
193 : !! seed_name= seed_name of files containing cg for all k-points to be used with MPI
194 : !! xred(3,natom)=reduced dimensionless atomic coordinates
195 : !!
196 : !! OUTPUT
197 : !! cm2: Inside sphere part of the overlap needed for constructing wannier function
198 : !!
199 : !! SIDE EFFECTS
200 : !! (only writing, printing)
201 : !!
202 : !! NOTES
203 : !! The mpi part will work with mlwfovlp but not for berryphase_new
204 : !!
205 : !! SOURCE
206 :
207 2944 : subroutine smatrix_pawinit(atindx1,cm2,cprj,ikpt1,ikpt2,isppol,&
208 2944 : & g1,gprimd,kpt,mband,mbandw,mkmem,mpi_enreg,&
209 2944 : & natom,nband,nkpt,nspinor,nsppol,ntypat,pawang,pawrad,pawtab,rprimd,&
210 2944 : & seed_name,typat,xred)
211 :
212 : !Arguments---------------------------
213 : !scalars
214 : integer,intent(in) :: ikpt1,ikpt2,isppol,mband,mbandw,mkmem,natom,nkpt,nspinor,nsppol
215 : integer,intent(in) :: ntypat
216 : character(len=fnlen) :: seed_name !seed names of files containing cg info used in case of MPI
217 : type(MPI_type),intent(in) :: mpi_enreg
218 : type(pawang_type),intent(in) :: pawang
219 :
220 : !arrays
221 : integer,intent(in) :: atindx1(natom),g1(3),nband(nsppol*nkpt),typat(natom)
222 : real(dp),intent(in) :: gprimd(3,3),kpt(3,nkpt),rprimd(3,3),xred(3,natom)
223 : real(dp),intent(inout) :: cm2(2,mbandw,mbandw)
224 : type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol)
225 : type(pawrad_type),intent(in) :: pawrad(ntypat)
226 : type(pawtab_type),intent(in) :: pawtab(ntypat)
227 :
228 : !Local variables---------------------------
229 : !scalars
230 : integer :: dummy
231 : integer :: iatom,iband1,iband2,icg1,icg2,idx1,idx2,ii
232 : integer :: ilmn,ios,iunit,ir
233 : integer :: iorder_cprj,isel,ispinor,itypat,j0lmn,jj,jlmn,klm,klmn,kln,ll,lm0,lmax
234 : integer :: lmin,lmn_size,max_lmn,mesh_size,mm,nband_k
235 : integer :: nprocs,spaceComm,rank !for mpi
236 : real(dp) :: arg,bnorm,delta,intg,ppi,ppr,qijbtemp,qijtot,x1
237 : real(dp) :: x2,xsum,xtemp,xx,yy,zz
238 : character(len=500) :: message
239 : character(len=fnlen) :: cprj_file !file containing cg info used in case of MPI
240 : logical::lfile
241 :
242 : !arrays
243 2944 : integer,allocatable :: dimcprj(:),nattyp_dum(:)
244 : real(dp),parameter :: ili(7)=(/zero,-one,zero,one,zero,-one,zero/)
245 : real(dp),parameter :: ilr(7)=(/one,zero,-one,zero,one,zero,-one/)
246 5888 : real(dp) :: bb(3),bb1(3),bbn(3),qijb(2),xcart(3,natom)
247 2944 : real(dp),allocatable :: ff(:),j_bessel(:,:),ylmb(:),ylmrgr_dum(:,:,:)
248 2944 : real(dp),allocatable :: sb_out(:)
249 : type(pawcprj_type),allocatable :: cprj_k1(:,:)
250 2944 : type(pawcprj_type),allocatable :: cprj_k2(:,:)
251 :
252 : ! *************************************************************************
253 :
254 : DBG_ENTER("COLL")
255 :
256 : !
257 : !Allocate cprj_k1 and cprj_k2
258 : !
259 8832 : ABI_MALLOC(dimcprj,(natom))
260 2944 : call pawcprj_getdim(dimcprj,natom,nattyp_dum,ntypat,typat,pawtab,'R')
261 :
262 2944 : nband_k=nband(ikpt1)
263 149632 : ABI_MALLOC(cprj_k1,(natom,nband_k*nspinor))
264 2944 : call pawcprj_alloc(cprj_k1,0,dimcprj)
265 :
266 2944 : nband_k=nband(ikpt2)
267 149632 : ABI_MALLOC(cprj_k2,(natom,nband_k*nspinor))
268 2944 : call pawcprj_alloc(cprj_k2,0,dimcprj)
269 2944 : ABI_FREE(dimcprj)
270 :
271 : !mpi initialization
272 2944 : spaceComm=MPI_enreg%comm_cell
273 2944 : nprocs=xmpi_comm_size(spaceComm)
274 2944 : rank=MPI_enreg%me_kpt
275 :
276 2944 : lfile=.false.
277 : !
278 : !write(std_out,*) "compute PAW overlap for k-points",ikpt1,ikpt2
279 6528 : do iatom=1,natom
280 : xcart(:,iatom)=rprimd(:,1)*xred(1,iatom)+&
281 : & rprimd(:,2)*xred(2,iatom)+&
282 17280 : & rprimd(:,3)*xred(3,iatom)
283 : end do
284 :
285 : !
286 : !Calculate indices icg1 and icg2
287 : !
288 2944 : icg1=0
289 6656 : do ii=1,isppol
290 3712 : ll=nkpt
291 3712 : if(ii==isppol) ll=ikpt1-1
292 144960 : do jj=1,ll
293 : ! MPI: cycle over kpts not treated by this node
294 138304 : if ( ABS(MPI_enreg%proc_distrb(jj,1,ii)-rank)/=0) CYCLE
295 : ! write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') jj,rank
296 142016 : icg1=icg1+nspinor*nband(jj+(ii-1)*nkpt)
297 : end do
298 : end do
299 2944 : icg2=0
300 6656 : do ii=1,isppol
301 3712 : ll=nkpt
302 3712 : if(isppol==ii) ll=ikpt2-1
303 144960 : do jj=1,ll
304 : ! MPI: cycle over kpts not treated by this node
305 138304 : if (ABS(MPI_enreg%proc_distrb(jj,1,ii)-rank)/=0) CYCLE
306 : ! write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') jj,rank
307 142016 : icg2=icg2+nspinor*nband(jj+(ii-1)*nkpt)
308 : end do
309 : end do
310 : !
311 : !MPI: if ikpt2 not found in this processor then
312 : !read info from an unformatted file
313 : !
314 2944 : if (nprocs>1) then
315 0 : if (ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-rank)/=0) then
316 0 : lfile=.true.
317 : !
318 : ! get maximum of lmn_size
319 : max_lmn=0
320 : do itypat=1,ntypat
321 : lmn_size=pawtab(itypat)%lmn_size
322 : if(lmn_size>max_lmn) max_lmn=lmn_size
323 : end do
324 : !
325 : ! get file name and open it
326 : !
327 0 : write(cprj_file,'(a,I5.5,".",I1)') trim(seed_name),ikpt2,isppol
328 0 : iunit=1000
329 : ! write(std_out,*)'reading file',trim(cprj_file)
330 0 : open (unit=iunit, file=cprj_file,form='unformatted',status='old',iostat=ios)
331 0 : if(ios /= 0) then
332 0 : write(message,*) " smatrix_pawinit: file",trim(cprj_file), "not found"
333 0 : ABI_ERROR(message)
334 : end if
335 : !
336 : ! start reading
337 0 : do ii=1,mband*nspinor
338 0 : do iatom=1,natom
339 0 : itypat=typat(iatom)
340 0 : lmn_size=pawtab(itypat)%lmn_size
341 0 : do ilmn=1,lmn_size
342 0 : read(iunit)(cprj_k2(iatom,ii)%cp(jj,ilmn),jj=1,2)
343 : end do !ilmn
344 : end do
345 : end do
346 : !
347 : ! close file
348 : !
349 0 : close (unit=iunit,iostat=ios)
350 0 : if(ios /= 0) then
351 0 : write(message,*) " smatrix_pawinit: error closing file ",trim(cprj_file)
352 0 : ABI_ERROR(message)
353 : end if
354 : !
355 : end if
356 : end if !mpi
357 :
358 : !Extract cprj_k1 and cprj_k2
359 : !these contain the projectors cprj for just one k-point (ikpt1 or ikpt2)
360 :
361 : !Extract cprj for k-point 1
362 2944 : iorder_cprj=0 !do not change the ordering of cprj
363 2944 : nband_k=nband(ikpt1)
364 2944 : dummy=1000 !index of file not implemented here, mkmem==0 not implemented
365 : call pawcprj_get(atindx1,cprj_k1,cprj,natom,1,icg1,ikpt1,iorder_cprj,isppol,&
366 : & mband,mkmem,natom,nband_k,nband_k,nspinor,nsppol,dummy,&
367 2944 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
368 :
369 : !Extract cprj for k-point 2
370 2944 : if( lfile .eqv. .false. ) then !if it was not already read above
371 : iorder_cprj=0 !do not change the ordering of cprj
372 2944 : nband_k=nband(ikpt2)
373 : dummy=1000 !index of file not implemented here, mkmem==0 not implemented
374 : call pawcprj_get(atindx1,cprj_k2,cprj,natom,1,icg2,ikpt2,iorder_cprj,isppol,&
375 : & mband,mkmem,natom,nband_k,nband_k,nspinor,nsppol,dummy,&
376 2944 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
377 : end if
378 :
379 : !DEBUG
380 : !if(ikpt2==2) then
381 : !if(ikpt1==1) then
382 : !do iband1=1,mbandw
383 : !do iatom=1,natom
384 : !itypat=typat(atindx1(iatom))
385 : !lmn_size=pawtab(itypat)%lmn_size
386 : !do ilmn=1,1!lmn_size
387 : !!write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',iatom,' ilmn ',ilmn,' cprj',cprj(iatom,iband1+icg1)%cp(:,ilmn)
388 : !write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',atindx1(iatom),' ilmn ',ilmn,' cprj',cprj(atindx1(iatom),iband1+icg1)%cp(:,ilmn)
389 : !write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',iatom,' ilmn ',ilmn,' cprj_k1 ',cprj_k1(iatom,iband1)%cp(:,ilmn)
390 : !end do
391 : !end do
392 : !end do
393 : !end if
394 : !end if
395 : !NDDEBUG
396 :
397 : !!!!!!!!!!!!!!!!!
398 : !--- Compute intermediate quantities: "b" vector=k2-k1 and its
399 : !normalized value: bbn (and its norm: bnorm)
400 : !compute also Ylm(b).
401 8832 : ABI_MALLOC(ylmb,(pawang%l_size_max*pawang%l_size_max))
402 2944 : ABI_MALLOC(ylmrgr_dum,(1,1,0))
403 11776 : bb(:)=kpt(:,ikpt2)-kpt(:,ikpt1)+g1(:)
404 : bb1=bb
405 2944 : xx=gprimd(1,1)*bb(1)+gprimd(1,2)*bb(2)+gprimd(1,3)*bb(3)
406 2944 : yy=gprimd(2,1)*bb(1)+gprimd(2,2)*bb(2)+gprimd(2,3)*bb(3)
407 2944 : zz=gprimd(3,1)*bb(1)+gprimd(3,2)*bb(2)+gprimd(3,3)*bb(3)
408 2944 : bnorm=two_pi*dsqrt(xx**2+yy**2+zz**2)
409 2944 : if(bnorm<tol8) then
410 : ! write(std_out,*) "WARNING: bnorm=",bnorm
411 0 : bbn(:)=zero
412 : else
413 2944 : xx=xx*two_pi
414 2944 : yy=yy*two_pi
415 2944 : zz=zz*two_pi
416 : bb(1)=xx
417 : bb(2)=yy
418 : bb(3)=zz
419 2944 : bbn(1)=xx/bnorm
420 2944 : bbn(2)=yy/bnorm
421 2944 : bbn(3)=zz/bnorm
422 : end if
423 :
424 : !debug bbn=0
425 : !debug bnorm=0
426 : !bbn has to ne normalized
427 2944 : call initylmr(pawang%l_size_max,0,1,(/one/),1,bbn(:),ylmb(:),ylmrgr_dum)
428 : !write(std_out,*) "ylmb(:)",ylmb(:)
429 : !write(std_out,*) pawang%l_size_max
430 : !write(std_out,*) "bbn",bbn(:)
431 : !write(std_out,*) "xx,yy,zz",xx,yy,zz
432 : !write(std_out,*) "bnorm",bnorm
433 2944 : ABI_FREE(ylmrgr_dum)
434 :
435 : !------- First Compute Qij(b)-
436 8832 : ABI_MALLOC(sb_out, (pawang%l_size_max))
437 2473472 : cm2=zero
438 :
439 6528 : do iatom=1,natom
440 3584 : itypat=typat(iatom)
441 3584 : lmn_size=pawtab(itypat)%lmn_size
442 : ! --- en coordonnnes reelles cartesiennes (espace reel)
443 : ! --- first radial part(see pawinit)
444 3584 : mesh_size=pawtab(itypat)%mesh_size
445 14336 : ABI_MALLOC(j_bessel,(mesh_size,pawang%l_size_max))
446 :
447 :
448 : ! --- compute bessel function for (br) for all angular momenta necessary
449 : ! --- and for all value of r.
450 : ! --- they are needed for radial part
451 : ! --- of the integration => j_bessel(ir,:)
452 4153600 : do ir=1,mesh_size
453 4150016 : arg=bnorm*pawrad(itypat)%rad(ir)
454 4150016 : call sbf8(pawang%l_size_max,arg,sb_out)
455 21701120 : j_bessel(ir,:) = sb_out
456 : end do
457 :
458 : ! do jlmn=1,pawang%l_size_max
459 : ! write(665,*) "j_bessel",j_bessel(1:mesh_size,jlmn)
460 : ! enddo
461 : ! write(std_out,*) "bessel function computed"
462 : ! --- Compute \Sum b.R=xsum for future use
463 : xtemp=zero
464 : do mm=1,3
465 : xtemp=xtemp+xred(mm,iatom)*bb1(mm)
466 : end do
467 : xtemp=xtemp*two_pi
468 : xsum=zero
469 14336 : do mm=1,3
470 14336 : xsum=xsum+xcart(mm,iatom)*bbn(mm)*bnorm
471 : end do
472 : ! write(std_out,*)'xsum',xsum,xtemp,lmn_size
473 :
474 : ! --- Loop on jlmn and ilmn
475 : qijtot=zero
476 56576 : do jlmn=1,lmn_size
477 52992 : j0lmn=jlmn*(jlmn-1)/2
478 513920 : do ilmn=1,jlmn
479 :
480 457344 : klmn=j0lmn+ilmn
481 457344 : klm=pawtab(itypat)%indklmn(1,klmn);kln=pawtab(itypat)%indklmn(2,klmn)
482 457344 : lmin=pawtab(itypat)%indklmn(3,klmn);lmax=pawtab(itypat)%indklmn(4,klmn)
483 : ! --- Sum over angular momenta
484 : ! --- compute radial part integration for each angular momentum => intg
485 : ! --- (3j) symbols follows the rule: l belongs to abs(li-lj), li+lj.
486 457344 : qijb=zero
487 457344 : do ll=lmin,lmax,2
488 946048 : lm0=ll*ll+ll+1
489 2838144 : ABI_MALLOC(ff,(mesh_size))
490 : ff(1:mesh_size)=(pawtab(itypat)%phiphj(1:mesh_size,kln)&
491 : & -pawtab(itypat)%tphitphj(1:mesh_size,kln))&
492 1014861440 : & *j_bessel(1:mesh_size,ll+1)
493 946048 : call simp_gen(intg,ff,pawrad(itypat))
494 946048 : ABI_FREE(ff)
495 946048 : qijbtemp=zero
496 5246208 : do mm=-ll,ll
497 4300160 : isel=pawang%gntselect(lm0+mm,klm)
498 4300160 : if (isel>0) qijbtemp=qijbtemp&
499 1865728 : & +pawang%realgnt(isel)*ylmb(lm0+mm)
500 : end do ! mm
501 : ! --- compute angular part with a summation
502 : ! --- qijb =\sum_{lm} intg(lm)*qijbtemp
503 946048 : qijb(1)=qijb(1) +intg*qijbtemp*ilr(ll+1)
504 1892096 : qijb(2)=qijb(2) +intg*qijbtemp*ili(ll+1)
505 : ! if(ilmn==jlmn) write(std_out,*) "intg, qij",intg,qijbtemp
506 : end do ! ll
507 :
508 : ! --- Add exp(-i.b*R) for each atom.
509 : if(ilmn==jlmn) qijtot=qijtot+qijb(1)
510 : ! if(ilmn==jlmn) write(std_out,*) "qijtot",qijtot
511 457344 : x1=qijb(1)*dcos(-xsum)-qijb(2)*dsin(-xsum)
512 457344 : x2=qijb(1)*dsin(-xsum)+qijb(2)*dcos(-xsum)
513 : ! x1 x2 necessary to avoid changing qijb(1) before
514 : ! computing qijb(2)
515 457344 : qijb(1)=x1
516 457344 : qijb(2)=x2 !
517 : ! if(ilmn==jlmn) write(std_out,*) "qij",jlmn,ilmn,qijb(1),qijb(2)
518 :
519 7917696 : do iband1=1,mbandw ! limite inferieure a preciser
520 143427456 : do iband2=1,mbandw
521 : ppr=0.d0
522 : ppi=0.d0
523 346770432 : do ispinor=1,nspinor
524 211207680 : idx1=iband1*nspinor-(nspinor-ispinor)
525 211207680 : idx2=iband2*nspinor-(nspinor-ispinor) !to take into account spinors
526 : ! write(std_out,*) "iband2",iband2
527 : ! product of (a1+ia2)*(b1-ib2) (minus sign because conjugated)
528 : ppr=ppr+&
529 : ! real part a_1*b_1+a_2*b_2
530 : & cprj_k1(iatom,idx1)%cp(1,ilmn)*cprj_k2(iatom,idx2)%cp(1,jlmn)+&
531 : & cprj_k1(iatom,idx1)%cp(2,ilmn)*cprj_k2(iatom,idx2)%cp(2,jlmn)+&
532 : ! & cprj(iatom,idx1+icg1)%cp(1,ilmn)*cprj(iatom,idx2+icg2)%cp(1,jlmn)+&
533 : ! & cprj(iatom,idx1+icg1)%cp(2,ilmn)*cprj(iatom,idx2+icg2)%cp(2,jlmn)+&
534 : ! add term on the other triangle of the matrix
535 : ! qij is the same for this part because phi are real.
536 : & cprj_k1(iatom,idx1)%cp(1,jlmn)*cprj_k2(iatom,idx2)%cp(1,ilmn)+&
537 211207680 : & cprj_k1(iatom,idx1)%cp(2,jlmn)*cprj_k2(iatom,idx2)%cp(2,ilmn)
538 : ! & cprj(iatom,idx1+icg1)%cp(1,jlmn)*cprj(iatom,idx2+icg2)%cp(1,ilmn)+&
539 : ! & cprj(iatom,idx1+icg1)%cp(2,jlmn)*cprj(iatom,idx2+icg2)%cp(2,ilmn)
540 : ppi=ppi+&
541 : ! imaginary part a_1*b_2-a_2*b_1
542 : & cprj_k1(iatom,idx1)%cp(1,ilmn)*cprj_k2(iatom,idx2)%cp(2,jlmn)-&
543 : & cprj_k1(iatom,idx1)%cp(2,ilmn)*cprj_k2(iatom,idx2)%cp(1,jlmn)+&
544 : ! & cprj(iatom,idx1+icg1)%cp(1,ilmn)*cprj(iatom,idx2+icg2)%cp(2,jlmn)-&
545 : ! & cprj(iatom,idx1+icg1)%cp(2,ilmn)*cprj(iatom,idx2+icg2)%cp(1,jlmn)+&
546 : ! add term on the other triangle of the matrix
547 : & cprj_k1(iatom,idx1)%cp(1,jlmn)*cprj_k2(iatom,idx2)%cp(2,ilmn)-&
548 346770432 : & cprj_k1(iatom,idx1)%cp(2,jlmn)*cprj_k2(iatom,idx2)%cp(1,ilmn)
549 : ! & cprj(iatom,idx1+icg1)%cp(1,jlmn)*cprj(iatom,idx2+icg2)%cp(2,ilmn)-&
550 : ! & cprj(iatom,idx1+icg1)%cp(2,jlmn)*cprj(iatom,idx2+icg2)%cp(1,ilmn)
551 : end do !ispinor
552 : !
553 : ! delta: diagonal terms are counted twice ! so
554 : ! we need a 0.5 factor for diagonal elements.
555 135562752 : delta=one
556 : ! write(std_out,*) "ppr and ppi computed",ikpt1,ikpt2,iband1,iband2
557 135562752 : if(ilmn==jlmn) delta=half
558 : cm2(1,iband1,iband2)= cm2(1,iband1,iband2)+ &
559 135562752 : & (qijb(1)*ppr-qijb(2)*ppi)*delta
560 : cm2(2,iband1,iband2)= cm2(2,iband1,iband2)+ &
561 142970112 : & (qijb(2)*ppr+qijb(1)*ppi)*delta
562 : end do ! iband2
563 : end do ! iband1
564 :
565 : end do ! ilmn
566 : end do ! jlmn
567 : ! write(std_out,*) "final qijtot",qijtot
568 6528 : ABI_FREE(j_bessel)
569 : end do ! iatom
570 :
571 2944 : ABI_FREE(sb_out)
572 2944 : ABI_FREE(ylmb)
573 2944 : call pawcprj_free(cprj_k1)
574 2944 : call pawcprj_free(cprj_k2)
575 76416 : ABI_FREE(cprj_k1)
576 76416 : ABI_FREE(cprj_k2)
577 :
578 : DBG_EXIT("COLL")
579 :
580 2944 : end subroutine smatrix_pawinit
581 : !!***
582 :
583 : !----------------------------------------------------------------------
584 :
585 : !!****f* m_paw_overlap/smatrix_k_paw
586 : !! NAME
587 : !! smatrix_k_paw
588 : !!
589 : !! FUNCTION
590 : !!
591 : !! INPUTS
592 : !! cprj_k (pawcprj_type) :: cprj for occupied bands at point k
593 : !! cprj_kb :: cprj for occupied bands at point k+b
594 : !! dtefield :: structure referring to all efield and berry's phase variables
595 : !! kdir :: integer giving direction along which overlap is computed for ket
596 : !! kfor :: integer indicating whether to compute forward (1) or backward (2)
597 : !! along kpt string
598 : !! natom :: number of atoms in cell
599 : !! typat :: typat(natom) type of each atom
600 : !!
601 : !! OUTPUT
602 : !! smat_k_paw :: array of the on-site PAW parts of the overlaps between Bloch states at points
603 : !! k and k+b, for the various pairs of bands, that is, the on-site part of
604 : !! <u_nk|u_mk+b>
605 : !!
606 : !! SIDE EFFECTS
607 : !!
608 : !! NOTES
609 : !! This routine assumes that the cprj are not explicitly ordered by
610 : !! atom type.
611 : !!
612 : !! SOURCE
613 :
614 1139688 : subroutine smatrix_k_paw(cprj_k,cprj_kb,dtefield,kdir,kfor,mband,natom,smat_k_paw,typat)
615 :
616 : !Arguments---------------------------
617 : !scalars
618 : integer,intent(in) :: kdir,kfor,mband,natom
619 : type(efield_type),intent(in) :: dtefield
620 : type(pawcprj_type),intent(in) :: cprj_k(natom,dtefield%nspinor*mband)
621 : type(pawcprj_type),intent(in) :: cprj_kb(natom,dtefield%nspinor*mband)
622 :
623 : !arrays
624 : integer,intent(in) :: typat(natom)
625 : real(dp),intent(out) :: smat_k_paw(2,dtefield%mband_occ,dtefield%mband_occ)
626 :
627 : !Local variables---------------------------
628 : !scalars
629 : integer :: iatom,iband,ibs,ilmn,ispinor,itypat
630 : integer :: jband,jbs,jlmn,klmn,nspinor
631 : complex(dp) :: cpk,cpkb,cterm,paw_onsite
632 : ! *************************************************************************
633 :
634 : !initialize smat_k_paw
635 60403464 : smat_k_paw(:,:,:) = zero
636 :
637 : nspinor = dtefield%nspinor
638 :
639 3419064 : do iatom = 1, natom
640 2279376 : itypat = typat(iatom)
641 :
642 21654072 : do ilmn=1,dtefield%lmn_size(itypat)
643 166394448 : do jlmn=1,dtefield%lmn_size(itypat)
644 145880064 : klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn)
645 : paw_onsite = cmplx(dtefield%qijb_kk(1,klmn,iatom,kdir),&
646 145880064 : & dtefield%qijb_kk(2,klmn,iatom,kdir))
647 145880064 : if (kfor > 1) paw_onsite = conjg(paw_onsite)
648 309995136 : do ispinor = 1, nspinor
649 875280384 : do iband = 1, dtefield%mband_occ
650 3063481344 : do jband = 1, dtefield%mband_occ
651 2334081024 : ibs = nspinor*(iband-1) + ispinor
652 2334081024 : jbs = nspinor*(jband-1) + ispinor
653 2334081024 : cpk=cmplx(cprj_k(iatom,ibs)%cp(1,ilmn),cprj_k(iatom,ibs)%cp(2,ilmn))
654 2334081024 : cpkb=cmplx(cprj_kb(iatom,jbs)%cp(1,jlmn),cprj_kb(iatom,jbs)%cp(2,jlmn))
655 2334081024 : cterm = conjg(cpk)*paw_onsite*cpkb
656 2334081024 : smat_k_paw(1,iband,jband) = smat_k_paw(1,iband,jband)+dreal(cterm)
657 2917601280 : smat_k_paw(2,iband,jband) = smat_k_paw(2,iband,jband)+dimag(cterm)
658 : end do ! end loop over jband
659 : end do ! end loop over iband
660 : end do ! end loop over ispinor
661 : end do ! end loop over ilmn
662 : end do ! end loop over jlmn
663 :
664 : end do ! end loop over atoms
665 :
666 1139688 : end subroutine smatrix_k_paw
667 : !!***
668 :
669 : !----------------------------------------------------------------------
670 :
671 : !!****f* m_paw_overlap/qijb_kk
672 : !! NAME
673 : !! qijb_kk
674 : !!
675 : !! FUNCTION
676 : !! Routine which computes PAW onsite part of wavefunction overlap for Bloch
677 : !! functions at two k-points k and k+b. These
678 : !! quantities are used in PAW-based computations of polarization and magnetization.
679 : !!
680 : !! INPUTS
681 : !! dkvecs(3) :: $\Delta k$ input vector
682 : !! expibi(2,my_natom,3) :: phase factors at each atomic site for given k offset
683 : !! gprimd(3,3)=dimensioned primitive translations of reciprocal lattice
684 : !! lmn2max :: lmnmax*(lmnmax+1)/2
685 : !! natom=number of atoms in unit cell
686 : !! ntypat=number of types of atoms in unit cell
687 : !! pawang <type(pawang_type)>=paw angular mesh and related data
688 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
689 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
690 : !! typat=typat(natom) list of atom types
691 : !!
692 : !! OUTPUT
693 : !! calc_qijb(2,lmn2max,natom) :: PAW on-site overlaps of wavefunctions at neighboring
694 : !! k point
695 : !!
696 : !! SIDE EFFECTS
697 : !!
698 : !! NOTES
699 : !! this function computes the on-site data for the PAW version of
700 : !! <u_nk|u_mk+b>, that is, two Bloch vectors at two different k points.
701 : !!
702 : !! SOURCE
703 :
704 90 : subroutine qijb_kk(calc_qijb,dkvecs,expibi,gprimd,lmn2max,natom,ntypat,&
705 90 : & pawang,pawrad,pawtab,typat)
706 :
707 : !Arguments---------------------------
708 : !scalars
709 : integer,intent(in) :: lmn2max,natom,ntypat
710 : type(pawang_type),intent(in) :: pawang
711 : real(dp),intent(out) :: calc_qijb(2,lmn2max,natom)
712 : !arrays
713 : integer,intent(in) :: typat(natom)
714 : real(dp),intent(in) :: dkvecs(3),expibi(2,natom),gprimd(3,3)
715 : type(pawrad_type),intent(in) :: pawrad(ntypat)
716 : type(pawtab_type),intent(in) :: pawtab(ntypat)
717 :
718 : !Local variables---------------------------
719 : !scalars
720 : integer :: iatom,ir,isel,itypat
721 : integer :: klm,kln,klmn,lbess,lbesslm,lmin,lmax,mbess,mesh_size
722 : integer :: ylmr_normchoice,ylmr_npts,ylmr_option
723 : real(dp) :: arg,bessg,bnorm,intg,rterm
724 : complex(dp) :: cterm,etb,ifac
725 : !arrays
726 : real(dp) :: bb(3),bbn(3),bcart(3),ylmgr(1,1,0),ylmr_nrm(1)
727 90 : real(dp),allocatable :: ff(:),j_bessel(:,:),ylmb(:),sb_out(:)
728 : ! the following is (i)^L mod 4.
729 : complex(dp),dimension(0:3) :: il(0:3)=(/cone,j_dpc,-cone,-j_dpc/)
730 : ! *************************************************************************
731 :
732 19710 : calc_qijb(:,:,:) = zero
733 :
734 90 : ylmr_normchoice = 0 ! input to initylmr are normalized
735 90 : ylmr_npts = 1 ! only 1 point to compute in initylmr
736 90 : ylmr_nrm(1) = one ! weight of normed point for initylmr
737 90 : ylmr_option = 1 ! compute only ylm's in initylmr
738 :
739 270 : ABI_MALLOC(sb_out, (pawang%l_size_max))
740 :
741 270 : do iatom = 1, natom
742 :
743 180 : itypat = typat(iatom)
744 180 : mesh_size = pawtab(itypat)%mesh_size
745 :
746 720 : ABI_MALLOC(j_bessel,(mesh_size,pawang%l_size_max))
747 540 : ABI_MALLOC(ff,(mesh_size))
748 540 : ABI_MALLOC(ylmb,(pawang%l_size_max*pawang%l_size_max))
749 :
750 : ! here is exp(-i b.R) for current atom: recall storage in expibi
751 180 : etb = cmplx(expibi(1,iatom),expibi(2,iatom))
752 :
753 : ! note the definition used for the k-dependence of the PAW basis functions:
754 : !$|\phi_{i,k}\rangle = exp(-i k\cdot r)|\phi_i\rangle
755 : ! see Umari, Gonze, and Pasquarello, PRB 69,235102 [[cite:Umari2004]], Eq. 23. Thus the k-vector on the
756 : ! bra side enters as k, while on the ket side it enters as -k.
757 720 : bb(:) = -dkvecs(:)
758 :
759 : ! reference bb to cartesian axes
760 2340 : bcart(1:3)=MATMUL(gprimd(1:3,1:3),bb(1:3))
761 :
762 : ! bbn is b-hat (the unit vector in the b direction)
763 720 : bnorm=dsqrt(dot_product(bcart,bcart))
764 720 : bbn(:) = bcart(:)/bnorm
765 :
766 : ! as an argument to the bessel function, need 2pi*b*r = 1 so b is re-normed to two_pi
767 180 : bnorm = two_pi*bnorm
768 92733 : do ir=1,mesh_size
769 92553 : arg=bnorm*pawrad(itypat)%rad(ir)
770 92553 : call sbf8(pawang%l_size_max,arg,sb_out) ! spherical bessel functions at each mesh point
771 370392 : j_bessel(ir,:) = sb_out
772 : end do ! end loop over mesh
773 :
774 : ! compute Y_LM(b) here
775 180 : call initylmr(pawang%l_size_max,ylmr_normchoice,ylmr_npts,ylmr_nrm,ylmr_option,bbn,ylmb(:),ylmgr)
776 :
777 6660 : do klmn = 1, pawtab(itypat)%lmn2_size
778 6480 : klm =pawtab(itypat)%indklmn(1,klmn)
779 6480 : kln =pawtab(itypat)%indklmn(2,klmn)
780 6480 : lmin=pawtab(itypat)%indklmn(3,klmn)
781 6480 : lmax=pawtab(itypat)%indklmn(4,klmn)
782 16920 : do lbess = lmin, lmax, 2 ! only possible choices for L s.t. Gaunt integrals
783 : ! will be non-zero
784 10260 : ifac = il(mod(lbess,4))
785 46440 : do mbess = -lbess, lbess
786 29700 : lbesslm = lbess*lbess+lbess+mbess+1
787 29700 : isel=pawang%gntselect(lbesslm,klm)
788 39960 : if (isel > 0) then
789 9180 : bessg = pawang%realgnt(isel)
790 : ff(1:mesh_size)=(pawtab(itypat)%phiphj(1:mesh_size,kln)&
791 : & -pawtab(itypat)%tphitphj(1:mesh_size,kln))&
792 4729383 : & *j_bessel(1:mesh_size,lbess+1)
793 9180 : call simp_gen(intg,ff,pawrad(itypat))
794 9180 : rterm = four_pi*bessg*intg*ylmb(lbesslm)
795 9180 : cterm = etb*ifac*rterm
796 : calc_qijb(1,klmn,iatom) = &
797 9180 : & calc_qijb(1,klmn,iatom) + dreal(cterm)
798 : calc_qijb(2,klmn,iatom) = &
799 9180 : & calc_qijb(2,klmn,iatom) + dimag(cterm)
800 :
801 : end if ! end selection on non-zero Gaunt factors
802 : end do ! end loop on mbess = -lbess, lbess
803 : end do ! end loop on lmin-lmax bessel l values
804 : end do ! end loop on lmn2_size klmn basis pairs
805 :
806 180 : ABI_FREE(j_bessel)
807 180 : ABI_FREE(ff)
808 270 : ABI_FREE(ylmb)
809 : end do ! end loop over atoms
810 :
811 90 : ABI_FREE(sb_out)
812 :
813 90 : end subroutine qijb_kk
814 : !!***
815 :
816 : !----------------------------------------------------------------------
817 :
818 : !!****f* m_paw_overlap/expibi
819 : !! NAME
820 : !! expibi
821 : !!
822 : !! FUNCTION
823 : !! Routine that computes exp(i (-b_ket).R) at each site.
824 : !!
825 : !! INPUTS
826 : !! dkvecs(3) :: $\Delta k$ increment
827 : !! natom :: number of atoms in unit cell
828 : !! xred(natom,3) :: reduced coordinates of atoms in unit cell
829 : !!
830 : !! OUTPUT
831 : !! calc_expibi(2,natom) :: phase factors at each atom for vector shift
832 : !!
833 : !! SIDE EFFECTS
834 : !!
835 : !! NOTES
836 : !!
837 : !! SOURCE
838 :
839 90 : subroutine expibi(calc_expibi,dkvecs,natom,xred)
840 :
841 : !Arguments---------------------------
842 : !scalars
843 : integer,intent(in) :: natom
844 : real(dp),intent(out) :: calc_expibi(2,natom)
845 : !arrays
846 : real(dp),intent(in) :: dkvecs(3),xred(3,natom)
847 :
848 : !Local variables---------------------------
849 : !scalars
850 : integer :: iatom
851 : real(dp) :: bdotr
852 :
853 : ! *************************************************************************
854 :
855 630 : calc_expibi(:,:) = zero
856 :
857 : !calc_expibi(2,natom)
858 : !used for PAW field calculations (distributed over atomic sites)
859 : !stores the on-site phase factors arising from
860 : !$\langle\phi_{i,k}|\phi_{j,k+\sigma_k k_k}\rangle$
861 : !where $\sigma = \pm 1$. These overlaps arise in various Berry
862 : !phase calculations of electric and magnetic polarization. The on-site
863 : !phase factor is $\exp[-i\sigma_k k_k)\cdot I]$ where
864 : !$I$ is the nuclear position.
865 :
866 270 : do iatom = 1, natom
867 :
868 : ! note the definition used for the k-dependence of the PAW basis functions:
869 : !$|\phi_{i,k}\rangle = exp(-i k\cdot r)|\phi_i\rangle
870 : ! see Umari, Gonze, and Pasquarello, PRB 69,235102 [[cite:Umari2004]] Eq. 23.
871 720 : bdotr = DOT_PRODUCT(xred(1:3,iatom),-dkvecs(1:3))
872 : ! here is exp(i b.R) for the given site
873 180 : calc_expibi(1,iatom) = cos(two_pi*bdotr)
874 270 : calc_expibi(2,iatom) = sin(two_pi*bdotr)
875 :
876 : end do ! end loop on natom
877 :
878 90 : end subroutine expibi
879 : !!***
880 :
881 : !----------------------------------------------------------------------
882 :
883 : END MODULE m_paw_overlap
884 : !!***
|
