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