Branch data Line data Source code
1 : : !!****m* ABINIT/m_wfd_optic
2 : : !! NAME
3 : : !! m_wfd_optic
4 : : !!
5 : : !! FUNCTION
6 : : !! Functions to compute optical matrix elements using the wavefunction descriptor.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2008-2026 ABINIT group (MG)
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_wfd_optic
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : : use m_xmpi
28 : :
29 : : use defs_datatypes, only : pseudopotential_type
30 : : use m_hide_lapack, only : matrginv
31 : : use m_bz_mesh, only : kmesh_t
32 : : use m_crystal, only : crystal_t
33 : : use m_vkbr, only : vkbr_t, vkbr_free, vkbr_init, nc_ihr_comm
34 : : use m_ebands, only : ebands_t
35 : : use m_wfd, only : wfdgw_t, wave_t
36 : : use m_pawtab, only : pawtab_type
37 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
38 : : use m_paw_hr, only : pawhur_t, paw_ihr
39 : :
40 : : implicit none
41 : :
42 : : private
43 : : !!***
44 : :
45 : : public :: calc_optical_mels
46 : : !!***
47 : :
48 : : contains
49 : : !!***
50 : :
51 : : !!****f* ABINIT/calc_optical_mels
52 : : !! NAME
53 : : !! calc_optical_mels
54 : : !!
55 : : !! FUNCTION
56 : : !! Calculate all optical matrix elements in the BZ.
57 : : !!
58 : : !! INPUTS
59 : : !! lomo_spin(Wfd%nsppol)=Index of the lomo band for the different spins.
60 : : !! lomo_min,max_band=minimum and max band index to be calculated.
61 : : !! nkbz=Number of points in the full Brillouin zone.
62 : : !! inclvkb=if different from 0, [Vnl,r] is included in the calculation of the
63 : : !! matrix element of the velocity operator. No meaning for PAW (except for DFT+U)
64 : : !! qpt(3)
65 : : !! Kmesh<kmesh_t>=Info on the k-point sampling for wave functions.
66 : : !! Cryst<crystal_t>=Structure defining the crystalline structure.
67 : : !! KS_Bst<ebands_t>
68 : : !! Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data
69 : : !! Psps <pseudopotential_type>=variables related to pseudopotentials.
70 : : !! Hur(Cryst%natom*usepaw)<pawhur_t>=Only for PAW and DFT+U, quantities used to evaluate the commutator [H_u,r].
71 : : !! Wfd<wfdgw_t>=Handler for the wavefunctions.
72 : : !!
73 : : !! OUTPUT
74 : : !! opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,nsppol)=Matrix elements <c k|e^{+iqr}|v k>
75 : : !!
76 : : !! SOURCE
77 : :
78 : 180 : subroutine calc_optical_mels(Wfd,Kmesh,KS_Bst,Cryst,Psps,Pawtab,Hur,&
79 : 180 : & inclvkb,lomo_spin,lomo_min,max_band,nkbz,qpoint,opt_cvk)
80 : :
81 : : !Arguments ------------------------------------
82 : : !scalars
83 : : integer,intent(in) :: nkbz,inclvkb,lomo_min,max_band
84 : : type(kmesh_t),intent(in) :: Kmesh
85 : : type(crystal_t),intent(in) :: Cryst
86 : : type(pseudopotential_type),intent(in) :: Psps
87 : : type(ebands_t),intent(in) :: KS_Bst
88 : : type(wfdgw_t),target,intent(inout) :: Wfd
89 : : !arrays
90 : : integer,intent(in) :: lomo_spin(Wfd%nsppol)
91 : : real(dp),intent(in) :: qpoint(3)
92 : : complex(dp),intent(out) :: opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,Wfd%nsppol)
93 : : type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
94 : : type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw)
95 : :
96 : : !Local variables ------------------------------
97 : : !scalars
98 : : integer :: nsppol,usepaw,nspinor,comm,spin,npw_k,istwf_k,my_nbbp
99 : : integer :: ik_bz,ik_ibz,itim_k,isym_k,ib_c,ib_v,ierr,my_rank
100 : : real(dp) :: ediff
101 : : complex(dp) :: emcvk
102 : : character(len=500) :: msg
103 : 180 : type(vkbr_t) :: vkbr
104 : : type(wave_t),pointer :: wave_v, wave_c
105 : : !arrays
106 : 180 : integer,allocatable :: bbp_distrb(:,:)
107 : 180 : integer,contiguous, pointer :: kg_k(:,:)
108 : : real(dp) :: mat_dp(3,3),qrot(3),b1(3),b2(3),b3(3),kbz(3)
109 : 180 : complex(dp),allocatable :: ir_kibz(:,:,:,:,:)
110 : 180 : complex(gwp), contiguous, pointer :: ug_c(:),ug_v(:)
111 : 360 : complex(gwp) :: ihrc(3,Wfd%nspinor**2)
112 : 360 : logical :: bbp_mask(Wfd%mband,Wfd%mband)
113 : 180 : type(pawcprj_type),allocatable :: Cp_v(:,:),Cp_c(:,:)
114 : : !************************************************************************
115 : :
116 : 180 : call wrtout(std_out," Calculating optical matrix elements in the IBZ","COLL")
117 [ - + ]: 180 : ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded")
118 : :
119 : 180 : comm = Wfd%comm
120 : 180 : my_rank = Wfd%my_rank
121 : :
122 : 180 : nsppol = Wfd%nsppol
123 : 180 : nspinor = Wfd%nspinor
124 : 180 : usepaw = Wfd%usepaw
125 : :
126 [ + + ]: 180 : if (usepaw==1) then
127 [ - + + - : 84 : ABI_MALLOC(Cp_v,(Wfd%natom,nspinor))
- + + - +
- + - - +
- + + + +
+ ]
128 : 12 : call pawcprj_alloc(Cp_v,0,Wfd%nlmn_atm)
129 [ - + + - : 84 : ABI_MALLOC(Cp_c,(Wfd%natom,nspinor))
- + + - +
- + - - +
- + + + +
+ ]
130 : 12 : call pawcprj_alloc(Cp_c,0,Wfd%nlmn_atm)
131 : : end if
132 : :
133 [ - + - - ]: 180 : if (inclvkb==1.and.usepaw==0) then
134 : 0 : ABI_ERROR("inclvkb==1 not coded,using inclvkb==2")
135 : : end if
136 : : !
137 : : ! Calculate the matrix elements of ir in the IBZ.
138 [ - + + - : 1080 : ABI_MALLOC(ir_kibz,(3,lomo_min:max_band,lomo_min:max_band,Wfd%nkibz,nsppol))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
139 [ + + + + : 1001316 : ir_kibz=czero
+ + + + +
+ ]
140 : :
141 [ - + + - : 720 : ABI_MALLOC(bbp_distrb, (Wfd%mband,Wfd%mband))
- + + - +
- + - - +
- + ]
142 : :
143 [ + + ]: 366 : do spin=1,nsppol
144 [ + + ]: 7122 : do ik_ibz=1,Wfd%nkibz
145 : : !
146 : : ! Distribute the (b,b') entries.
147 [ + + + + : 660552 : bbp_mask=.FALSE.; bbp_mask(lomo_spin(spin):max_band,lomo_spin(spin):max_band)=.TRUE.
+ + + + ]
148 : 6756 : call wfd%distribute_bbp(ik_ibz,spin,"All",my_nbbp,bbp_distrb,bbp_mask=bbp_mask)
149 [ + - + + : 72762 : if (ALL(bbp_distrb/=my_rank)) CYCLE
+ + - + ]
150 : :
151 : 6756 : istwf_k = Wfd%istwfk(ik_ibz)
152 [ - + ]: 6756 : ABI_CHECK(istwf_k==1,"istwf_k/=1 not coded") ! KB stuff is missing.
153 : 6756 : npw_k = Wfd%npwarr(ik_ibz)
154 : 6756 : kg_k => Wfd%Kdata(ik_ibz)%kg_k
155 : :
156 [ + - + + ]: 6756 : if (inclvkb/=0.and.usepaw==0) then
157 : : ! Prepare term i <n,k|[Vnl,r]|n"k>
158 : 6714 : call vkbr_init(vkbr,Cryst,Psps,inclvkb,istwf_k,npw_k,Kmesh%ibz(:,ik_ibz),kg_k)
159 : : end if
160 : :
161 : : ! Note: spinorial case is not coded therefore we work with ihrc(:,1).
162 : : ! TODO: The lower triangle can be Reconstructed by symmetry.
163 [ + + ]: 46350 : do ib_v=lomo_spin(spin),max_band ! Loop over bands
164 [ + - + + : 94164 : if ( ALL(bbp_distrb(ib_v,:)/=my_rank) ) CYCLE
+ - ]
165 : :
166 [ - + ]: 39594 : ABI_CHECK(wfd%get_wave_ptr(ib_v, ik_ibz, spin, wave_v, msg) == 0, msg)
167 : 39594 : ug_v => wave_v%ug
168 [ + + ]: 39594 : if (usepaw==1) call wfd%get_cprj(ib_v,ik_ibz,spin,Cryst,Cp_v,sorted=.FALSE.)
169 : :
170 [ + + ]: 285000 : do ib_c=lomo_spin(spin),max_band
171 [ + + ]: 238650 : if (bbp_distrb(ib_v,ib_c)/=my_rank) CYCLE
172 [ - + ]: 200250 : ABI_CHECK(wfd%get_wave_ptr(ib_c, ik_ibz, spin, wave_c, msg) == 0, msg)
173 : 200250 : ug_c => wave_c%ug
174 : :
175 [ + + ]: 200250 : if (usepaw==0) then
176 : : ! Calculate matrix elements of i[H,r] for NC pseudopotentials.
177 : 198426 : ihrc = nc_ihr_comm(vkbr,cryst,psps,npw_k,nspinor,istwf_k,inclvkb,Kmesh%ibz(:,ik_ibz),ug_c,ug_v,kg_k)
178 : :
179 : : else
180 : : ! Matrix elements of i[H,r] for PAW.
181 : 1824 : call wfd%get_cprj(ib_c,ik_ibz,spin,Cryst,Cp_c,sorted=.FALSE.)
182 : :
183 : 1824 : ihrc = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug_c,ug_v,kg_k,Cp_c,Cp_v,HUr)
184 : : end if
185 : : !
186 : : ! Save matrix elements of i*r in the IBZ
187 : 200250 : ediff = KS_Bst%eig(ib_c,ik_ibz,spin) - KS_BSt%eig(ib_v,ik_ibz,spin)
188 [ + + ]: 200250 : if (ABS(ediff)<tol16) ediff=tol6 ! Treat a possible degeneracy between v and c.
189 [ + + ]: 840594 : ir_kibz(:,ib_c,ib_v,ik_ibz,spin) = ihrc(:,1)/ediff
190 : :
191 : : end do !ib_c
192 : : end do !ib_v
193 : :
194 : 6942 : call vkbr_free(vkbr)
195 : : end do !spin
196 : : end do !ik_ibz
197 : :
198 : : ! Collect results on each node.
199 : 180 : call xmpi_sum(ir_kibz,comm,ierr)
200 : :
201 [ - + ]: 180 : ABI_FREE(bbp_distrb)
202 : :
203 [ + + ]: 180 : if (usepaw==1) then
204 : 12 : call pawcprj_free(Cp_v)
205 [ + - + + : 36 : ABI_FREE(Cp_v)
- + - + -
+ ]
206 : 12 : call pawcprj_free(Cp_c)
207 [ + - + + : 36 : ABI_FREE(Cp_c)
- + - + -
+ ]
208 : : end if
209 : : !
210 : : ! ======================================================
211 : : ! ==== Calculate Fcv(kBZ) in the full Brilouin zone ====
212 : : ! ======================================================
213 : : !
214 : : ! Symmetrization of the matrix elements of the position operator.
215 : : ! <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
216 : : ! where S is one of the symrec operations in reciprocal space, R is the
217 : : ! corresponding operation in real space, \tau being the associated fractional translations.
218 : : !
219 : : ! q.Mcv( Sk) = S^{-1}q. Mcv(k)
220 : : ! q.Mcv(-Sk) = -S^{-1}q. CONJG(Mcv(k)) if time-reversal is used.
221 : :
222 [ + + ]: 720 : b1=Cryst%gprimd(:,1)*two_pi
223 [ + + ]: 720 : b2=Cryst%gprimd(:,2)*two_pi
224 [ + + ]: 720 : b3=Cryst%gprimd(:,3)*two_pi
225 : :
226 [ + + + + : 384126 : opt_cvk = czero
+ + + + ]
227 [ + + ]: 366 : do spin=1,nsppol
228 [ + + ]: 8862 : do ik_bz=1,nkbz
229 : : !
230 : : ! Get ik_ibz, and symmetries index from ik_bz.
231 : 8496 : call Kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k)
232 : :
233 [ + + + + ]: 110448 : mat_dp = DBLE(Cryst%symrec(:,:,isym_k))
234 : 8496 : call matrginv(mat_dp,3,3) ! Invert
235 [ + + + + : 135936 : qrot = (3-2*itim_k) * MATMUL(mat_dp,qpoint)
+ + ]
236 : :
237 [ + + ]: 60426 : do ib_v=lomo_spin(spin),max_band ! Loops over the bands C and V start
238 [ + + ]: 383760 : do ib_c=lomo_spin(spin),max_band
239 : : !if (ib_c==ib_v) CYCLE
240 : 323520 : emcvk = pdtqrc(qrot,ir_kibz(:,ib_c,ib_v,ik_ibz,spin),b1,b2,b3)
241 [ + + ]: 323520 : if (itim_k==2) emcvk = CONJG(emcvk)
242 : 375264 : opt_cvk(ib_c,ib_v,ik_bz,spin) = emcvk
243 : : end do !ib_c
244 : : end do !ib_v
245 : :
246 : : end do !ik_bz
247 : : end do !spin
248 : :
249 [ - + ]: 180 : ABI_FREE(ir_kibz)
250 : :
251 [ - + - + : 540 : call xmpi_barrier(comm)
- + - - -
- - - - +
- - - - -
- ]
252 : :
253 : : contains
254 : : !!***
255 : :
256 : : !!****f* ABINIT/pdtqrc
257 : : !! NAME
258 : : !! pdtqrc
259 : : !!
260 : : !! FUNCTION
261 : : !! Calculate the dot product of a real vector with a complex vector, where each is in terms of b1-b3
262 : : !!
263 : : !! INPUTS
264 : : !!
265 : : !! OUTPUT
266 : : !!
267 : : !! SOURCE
268 : :
269 : 323520 : pure function pdtqrc(R,C,b1,b2,b3)
270 : :
271 : : !Arguments ------------------------------------
272 : : !arrays
273 : : real(dp),intent(in) :: R(3),b1(3),b2(3),b3(3)
274 : : complex(dp),intent(in) :: C(3)
275 : : complex(dp) :: pdtqrc
276 : :
277 : : !Local variables ------------------------------
278 : : !scalars
279 : : integer :: ii
280 : : !************************************************************************
281 : :
282 : 323520 : pdtqrc=czero
283 [ + + ]: 1294080 : do ii=1,3
284 : : pdtqrc = pdtqrc + (R(1)*b1(ii)+R(2)*b2(ii)+R(3)*b3(ii)) * &
285 : 1294080 : & (C(1)*b1(ii)+C(2)*b2(ii)+C(3)*b3(ii))
286 : : end do
287 : :
288 : 323520 : end function pdtqrc
289 : : !!***
290 : :
291 : : end subroutine calc_optical_mels
292 : : !!***
293 : :
294 : : end module m_wfd_optic
295 : : !!***
|
