Line data Source code
1 : !****m* ABINIT/m_oscillators
2 : !! NAME
3 : !! m_oscillators
4 : !!
5 : !! FUNCTION
6 : !! This module contains procedures to calculate the oscillator matrix elements used in the GW code.
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_oscillators
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_fft
28 :
29 : use m_gwdefs, only : czero_gw
30 : use m_fstrings, only : toupper, itoa, sjoin
31 : use m_geometry, only : spinrot_cmat
32 : use m_hide_blas, only : xcopy
33 : use m_gsphere, only : gsphere_t
34 :
35 : implicit none
36 :
37 : private
38 : !!***
39 :
40 : !----------------------------------------------------------------------
41 :
42 : public :: rho_tw_g ! Calculate rhotwg(G) = <wfn1| exp(-i(q+G).r) |wfn2>
43 : public :: calc_wfwfg ! Calculate the Fourier transform of the product u_{bk}^*(r).u_{b"k}(r) at an arbitrary k in the BZ.
44 : public :: sym_rhotwgq0 ! Symmetrize the oscillator matrix elements in the BZ in the special case of q = 0.
45 : !!***
46 :
47 : !----------------------------------------------------------------------
48 :
49 : CONTAINS
50 :
51 : !!****f* m_oscillators/rho_tw_g
52 : !! NAME
53 : !! rho_tw_g
54 : !!
55 : !! FUNCTION
56 : !! Calculate rhotwg(G) = <wfn1| exp(-i(q+G).r) |wfn2>
57 : !!
58 : !! INPUTS
59 : !! dim_rtwg=Define the size of the output array rhotwg
60 : !! === for nspinor==1 ===
61 : !! dim_rtwg=1
62 : !! === for nspinor==2 ===
63 : !! dim_rtwg=1 if the sum of the matrix elements is wanted.
64 : !! dim_rtwg=2 if <up|up>, <dwn|dwn> matrix elements are required
65 : !! map2sphere= 1 to retrieve Fourier components indexed according to igfftg0.
66 : !! 0 to retrieve Fourier components indexed according to the FFT box.
67 : !! NOTE: If map2sphere==0 npwvec must be equal to nr
68 : !! use_padfft= Only compatible with map2sphere 1.
69 : !! 1 if matrix elements are calculated via zero-padded FFT.
70 : !! 0 R-->G Transform in done on the full FFT box.
71 : !! igfftg0(npwvec*map2sphere)=index of G-G_o in the FFT array for each G in the sphere.
72 : !! i1=1 if kbz1 = Sk1, 2 if kbz1 = -Sk_1 (k_1 is in the IBZ)
73 : !! i2=1 if kbz2 = Sk2, 2 if kbz2 = -Sk_2 (k_2 is in the IBZ)
74 : !! ktabr1(nr),ktabr2(nr)= tables R^-1(r-t) for the two k-points
75 : !! ktabp1,ktabp2 = phase factors for non-simmorphic symmetries e^{-i 2\pi kbz.\tau}
76 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
77 : !! npwvec=number of plane waves (in the sphere if map2sphere==1, in the FFT box if map2sphere==1)
78 : !! nr=number of FFT grid points
79 : !! ndat=Number of wavefunctions to transform.
80 : !! nspinor=number of spinorial components.
81 : !! spinrot1(4),spinrot2(4)=components of the spinor rotation matrix. See getspinrot
82 : !! wfn1(nr*nspinor*ndat),wfn2(nr*nspinor*ndat)=the two wavefunctions (periodic part)
83 : !! [nhat12(2,nr,nspinor**2*ndat)]=Compensation charge in real space to be added to \Psi_1^*\Psi_2 -- Only for PAW.
84 : !!
85 : !! OUTPUT
86 : !! rhotwg(npwvec)=density of a pair of states, in reciprocal space
87 : !!
88 : !! SOURCE
89 :
90 5694935 : subroutine rho_tw_g(nspinor, npwvec, nr, ndat, ngfft, map2sphere, use_padfft, igfftg0, gbound, &
91 5694935 : wfn1, i1, ktabr1, ktabp1, spinrot1, &
92 5694935 : wfn2, i2, ktabr2, ktabp2, spinrot2, &
93 5694935 : dim_rtwg, rhotwg) !& nhat12)
94 :
95 : !Arguments ------------------------------------
96 : !scalars
97 : integer,intent(in) :: i1,i2,npwvec,nr,nspinor,dim_rtwg,map2sphere,use_padfft,ndat
98 : complex(dp),intent(in) :: ktabp1, ktabp2
99 : !arrays
100 : integer,intent(in) :: gbound(:,:) !gbound(2*mgfft+8,2)
101 : integer,intent(in) :: igfftg0(npwvec*map2sphere),ngfft(18)
102 : integer,intent(in) :: ktabr1(nr),ktabr2(nr)
103 : real(dp),intent(in) :: spinrot1(4),spinrot2(4)
104 : complex(gwp),intent(in) :: wfn1(nr*nspinor*ndat),wfn2(nr*nspinor*ndat)
105 : complex(gwp),intent(out) :: rhotwg(npwvec*dim_rtwg*ndat)
106 : ! real(dp),optional,intent(in) :: nhat12(2,nr,nspinor**2*ndat)
107 :
108 : !Local variables-------------------------------
109 : !scalars
110 : integer :: fftcache0 = 0, gpu_option_0 = 0
111 : integer :: ig,igfft,iab,spad1,spad2,spad0,nx,ny,nz,ldx,ldy,ldz,mgfft
112 : type(fftbox_plan3_t) :: plan
113 : !arrays
114 : integer :: spinor_pad(2,4)
115 5694935 : complex(gwp),allocatable :: u12prod(:),cwavef1(:),cwavef2(:),cwork(:)
116 : ! *************************************************************************
117 :
118 5643847 : SELECT CASE (nspinor)
119 : CASE (1)
120 : ! Collinear case.
121 : call ts_usug_kkp_bz(npwvec,nr,ndat,ngfft,map2sphere,use_padfft,igfftg0,gbound,&
122 : wfn1,i1,ktabr1,ktabp1,&
123 5643847 : wfn2,i2,ktabr2,ktabp2,rhotwg)
124 :
125 : CASE (2)
126 : ! Spinorial case.
127 51088 : ABI_CHECK(ndat==1,"ndat != 1 not coded")
128 153264 : ABI_MALLOC(cwavef1, (nr * nspinor * ndat))
129 102176 : ABI_MALLOC(cwavef2, (nr * nspinor * ndat))
130 102176 : ABI_MALLOC(cwork, (nr * nspinor * ndat))
131 :
132 51088 : call rotate_spinor(i1, ktabr1, ktabp1, spinrot1, nr, nspinor, ndat, wfn1, cwork, cwavef1)
133 51088 : call rotate_spinor(i2, ktabr2, ktabp2, spinrot2, nr, nspinor, ndat, wfn2, cwork, cwavef2)
134 :
135 51088 : ABI_FREE(cwork)
136 153264 : ABI_MALLOC(u12prod, (nr))
137 :
138 459792 : spinor_pad = reshape([0, 0, nr, nr, 0, nr, nr, 0], [2, 4])
139 17654704 : rhotwg = czero_gw
140 153264 : do iab=1,2
141 102176 : spad1 = spinor_pad(1,iab); spad2=spinor_pad(2,iab)
142 :
143 817612352 : u12prod = GWPC_CONJG(cwavef1(spad1+1:spad1+nr)) * cwavef2(spad2+1:spad2+nr)
144 : ! Add compensation charge.
145 : !if (PRESENT(nhat12)) u12prod = u12prod + CMPLX(nhat12(1,:,iab),nhat12(2,:,iab))
146 :
147 102176 : spad0 = (iab-1)*npwvec
148 51088 : SELECT CASE (map2sphere)
149 : CASE (0)
150 : ! Need results on the full FFT box thus cannot use zero-padded FFT.
151 0 : call plan%init(ndat, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0)
152 0 : call plan%execute(u12prod, -1, ndat)
153 0 : call plan%free()
154 0 : if (dim_rtwg == 1) then
155 0 : rhotwg(1:npwvec) = rhotwg(1:npwvec) + u12prod
156 : else
157 0 : rhotwg(spad0+1:spad0+npwvec) = u12prod
158 : end if
159 :
160 : CASE (1)
161 : ! Need results on the G-sphere. Call zero-padded FFT routines if required.
162 102176 : if (use_padfft == 1) then
163 408704 : nx = ngfft(1); ny = ngfft(2); nz = ngfft(3); mgfft = maxval(ngfft(1:3))
164 102176 : ldx = nx; ldy = ny; ldz = nz
165 102176 : call fftpad(u12prod, ngfft, nx, ny, nz, ldx, ldy, ldz, ndat, mgfft, -1, gbound)
166 : else
167 0 : call plan%init(ndat, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0)
168 0 : call plan%execute(u12prod, -1, ndat)
169 0 : call plan%free()
170 : end if
171 :
172 : ! Have to map FFT to G-sphere.
173 102176 : if (dim_rtwg == 1) then
174 354816 : do ig=1,npwvec
175 342144 : igfft = igfftg0(ig)
176 : ! G-G0 belong to the FFT mesh.
177 354816 : if (igfft /= 0) rhotwg(ig) = rhotwg(ig) + u12prod(igfft)
178 : end do
179 : else
180 17522048 : do ig=1,npwvec
181 17432544 : igfft = igfftg0(ig)
182 : ! G-G0 belong to the FFT mesh.
183 17522048 : if (igfft /= 0) rhotwg(ig+spad0) = u12prod(igfft)
184 : end do
185 : end if
186 :
187 : CASE DEFAULT
188 102176 : ABI_BUG("Wrong map2sphere")
189 : END SELECT
190 : end do !iab
191 :
192 51088 : ABI_FREE(u12prod)
193 51088 : ABI_FREE(cwavef1)
194 51088 : ABI_FREE(cwavef2)
195 :
196 : CASE DEFAULT
197 0 : ABI_BUG('Wrong nspinor')
198 : END SELECT
199 :
200 45559480 : end subroutine rho_tw_g
201 : !!***
202 :
203 : !----------------------------------------------------------------------
204 :
205 : !!****f* m_oscillators/ts_usug_kkp_bz
206 : !! NAME
207 : !! ts_usug_kkp_bz
208 : !!
209 : !! FUNCTION
210 : !! Calculate usug(G) = <u1|exp(-i(q+G).r)|u2> for ndat pair of wavefunctions
211 : !! TODO: The routine is thread-safe hence it can be called within an OMP parallel region.
212 : !!
213 : !! INPUTS
214 : !! map2sphere= 1 to retrieve Fourier components indexed according to igfftg0.
215 : !! 0 to retrieve Fourier components indexed according to the FFT box.
216 : !! NOTE: If map2sphere==0 npw must be equal to nr
217 : !! use_padfft= Only compatible with map2sphere 1.
218 : !! 1 if matrix elements are calculated via zero-padded FFT.
219 : !! 0 R-->G Transform in done on the full FFT box.
220 : !! igfftg0(npw*map2sphere)=index of G-G_o in the FFT array for each G in the sphere.
221 : !! time1=1 if kbz1 = Sk1, 2 if kbz1 = -Sk_1 (k_1 is in the IBZ)
222 : !! time2=1 if kbz2 = Sk2, 2 if kbz2 = -Sk_2 (k_2 is in the IBZ)
223 : !! ktabr1(nr),ktabr2(nr)= tables R^-1(r-t) for the two k-points
224 : !! ktabp1,ktabp2 = phase factors for non-simmorphic symmetries e^{-i 2\pi kbz.\tau}
225 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
226 : !! npw=number of plane waves (in the sphere if map2sphere==1, in the FFT box if map2sphere==1)
227 : !! nr=number of FFT grid points
228 : !! ndat=Number of wavefunctions to transform.
229 : !! u1(nr*ndat),u2(nr*ndat)=the two wavefunctions (periodic part)
230 : !!
231 : !! OUTPUT
232 : !! usug(npw*ndat)=density of a pair of states, in reciprocal space
233 : !!
234 : !! SOURCE
235 :
236 5643847 : subroutine ts_usug_kkp_bz(npw, nr, ndat, ngfft, map2sphere, use_padfft, igfftg0, gbound, &
237 5643847 : u1, time1, ktabr1, ktabp1, &
238 5643847 : u2, time2, ktabr2, ktabp2, usug) !& nhat12)
239 :
240 : !Arguments ------------------------------------
241 : !scalars
242 : integer,intent(in) :: time1,time2,npw,nr,map2sphere,use_padfft,ndat
243 : complex(dp),intent(in) :: ktabp1,ktabp2
244 : !arrays
245 : integer,intent(in) :: gbound(:,:) !gbound(2*mgfft+8,2)
246 : integer,intent(in) :: igfftg0(npw*map2sphere),ngfft(18)
247 : integer,intent(in) :: ktabr1(nr),ktabr2(nr)
248 : complex(gwp),intent(in) :: u1(nr*ndat),u2(nr*ndat)
249 : complex(gwp),intent(out) :: usug(npw*ndat)
250 :
251 : !Local variables-------------------------------
252 : !scalars
253 : integer :: fftcache0 = 0, gpu_option_0 = 0
254 : integer :: nx,ny,nz,ldx,ldy,ldz,mgfft
255 : type(fftbox_plan3_t) :: plan
256 : !arrays
257 5643847 : complex(gwp),allocatable :: u12prod(:)
258 : ! *************************************************************************
259 :
260 : ! Form rho-twiddle(r) = u_1^*(r,b1,kbz1) u_2(r,b2,kbz2), to account for symmetries:
261 : !
262 : ! u(r,b,kbz) = e^{-2i\pi kibz.(R^{-1}t} u (R{^-1}(r-t), b, kibz)
263 : ! = e^{+2i\pi kibz.(R^{-1}t} u*({R^-1}(r-t), b, kibz) for time-reversal symmetry.
264 : !
265 16931541 : ABI_MALLOC(u12prod,(nr*ndat))
266 5643847 : call usur_kkp_bz(nr,ndat,time1,ktabr1,ktabp1,u1,time2,ktabr2,ktabp2,u2,u12prod)
267 :
268 : ! Add compensation charge.
269 : !if (PRESENT(nhat12)) u12prod = u1prod + CMPLX(nhat12(1,:,1),nhat12(2,:,1))
270 :
271 0 : SELECT CASE (map2sphere)
272 : CASE (0)
273 : ! Need results on the full FFT box thus cannot use zero-padded FFT.
274 0 : call plan%init(ndat, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0)
275 0 : call plan%execute(u12prod, -1, ndat)
276 0 : call plan%free()
277 0 : call xcopy(nr*ndat,u12prod,1,usug,1)
278 :
279 : CASE (1)
280 : ! Need results on the G-sphere. Call zero-padded FFT routines if required.
281 5643847 : if (use_padfft==1) then
282 22526960 : nx = ngfft(1); ny = ngfft(2); nz = ngfft(3); mgfft = MAXVAL(ngfft(1:3))
283 5631740 : ldx=nx; ldy=ny; ldz=nz
284 5631740 : call fftpad(u12prod,ngfft,nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,-1,gbound)
285 : else
286 12107 : call plan%init(ndat, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0)
287 12107 : call plan%execute(u12prod, -1, ndat)
288 12107 : call plan%free()
289 : end if
290 :
291 : ! From the FFT to the G-sphere.
292 5643847 : call gw_box2gsph(nr,ndat,npw,igfftg0,u12prod,usug)
293 :
294 : CASE DEFAULT
295 5643847 : ABI_BUG("Wrong map2sphere")
296 : END SELECT
297 :
298 5643847 : ABI_FREE(u12prod)
299 :
300 45150776 : end subroutine ts_usug_kkp_bz
301 : !!***
302 :
303 : !----------------------------------------------------------------------
304 :
305 : !!****f* m_oscillators/usur_kkp_bz
306 : !! NAME
307 : !! usur_kkp_bz
308 : !!
309 : !! FUNCTION
310 : !! Calculate u1_kbz^*(r) u2_kbz(r) in real space from the symmetric images in the IBZ.
311 : !! Does not support spinor wavefunctions.
312 : !!
313 : !! INPUTS
314 : !! nr=number of FFT grid points
315 : !! ndat=Number of wavefunctions to transform.
316 : !! u1(nr*ndat),u2(nr*ndat)=the two wavefunctions in the IBZ (periodic part)
317 : !! time1=1 if kbz1 = Sk1, 2 if kbz1 = -Sk_1 (k_1 is in the IBZ)
318 : !! time2=1 if kbz2 = Sk2, 2 if kbz2 = -Sk_2 (k_2 is in the IBZ)
319 : !! ktabr1(nr),ktabr2(nr)= tables R^-1(r-t) for the two k-points
320 : !! ktabp1,ktabp2 = phase factors for non-simmorphic symmetries e^{-i 2\pi kbz.\tau}
321 : !!
322 : !! OUTPUT
323 : !! u12prod(nr*dat) = u1_kbz^*(r) u2_kbz(r) for the ndat pairs.
324 : !!
325 : !! SOURCE
326 :
327 5643847 : subroutine usur_kkp_bz(nr, ndat, time1, ktabr1, ktabp1, u1, time2, ktabr2, ktabp2, u2, u12prod)
328 :
329 : !Arguments ------------------------------------
330 : !scalars
331 : integer,intent(in) :: nr,ndat,time1,time2
332 : complex(dp),intent(in) :: ktabp1,ktabp2
333 : !arrays
334 : integer,intent(in) :: ktabr1(nr),ktabr2(nr)
335 : complex(gwp),intent(in) :: u1(nr*ndat),u2(nr*ndat)
336 : complex(gwp),intent(out) :: u12prod(nr*ndat)
337 :
338 : !Local variables-------------------------------
339 : !scalars
340 : integer :: ir,dat,padat
341 : complex(gwp) :: my_ktabp1,my_ktabp2
342 : !arrays
343 5643847 : complex(gwp),allocatable :: u1_bz(:),u2_bz(:)
344 : ! *************************************************************************
345 :
346 : ! Form rho-twiddle(r)=u_1^*(r,b1,kbz1) u_2(r,b2,kbz2), to account for symmetries:
347 : ! u(r,b,kbz)=e^{-2i\pi kibz.(R^{-1}t} u (R{^-1}(r-t),b,kibz)
348 : ! =e^{+2i\pi kibz.(R^{-1}t} u*({R^-1}(r-t),b,kibz) for time-reversal
349 : !
350 16931541 : ABI_MALLOC(u1_bz,(nr*ndat))
351 11287694 : ABI_MALLOC(u2_bz,(nr*ndat))
352 :
353 5643847 : my_ktabp1 = ktabp1
354 5643847 : my_ktabp2 = ktabp2
355 :
356 5643847 : if (ndat==1) then
357 33854934128 : do ir=1,nr
358 33854934128 : u1_bz(ir) = u1(ktabr1(ir))*my_ktabp1
359 : end do
360 33854934128 : do ir=1,nr
361 33854934128 : u2_bz(ir) = u2(ktabr2(ir))*my_ktabp2
362 : end do
363 : else
364 : !$OMP PARALLEL PRIVATE(padat)
365 : !$OMP DO
366 0 : do dat=1,ndat
367 0 : padat = (dat-1)*nr
368 0 : do ir=1,nr
369 0 : u1_bz(ir+padat) = u1(ktabr1(ir)+padat)*my_ktabp1
370 : end do
371 : end do
372 : !$OMP END DO NOWAIT
373 : !$OMP DO
374 0 : do dat=1,ndat
375 0 : padat = (dat-1)*nr
376 0 : do ir=1,nr
377 0 : u2_bz(ir+padat) = u2(ktabr2(ir)+padat)*my_ktabp2
378 : end do
379 : end do
380 : !$OMP END DO NOWAIT
381 : !$OMP END PARALLEL
382 : end if
383 :
384 : ! Treat time-reversal.
385 10920594 : SELECT CASE (time1)
386 : CASE (1)
387 5276747 : if (ndat==1) then
388 5276747 : if (time2==1) then
389 32558152608 : do ir=1,nr
390 32558152608 : u12prod(ir) = GWPC_CONJG(u1_bz(ir)) * u2_bz(ir)
391 : end do
392 70980 : else if (time2==2) then
393 191450328 : do ir=1,nr
394 191450328 : u12prod(ir) = GWPC_CONJG(u1_bz(ir)) * GWPC_CONJG(u2_bz(ir))
395 : end do
396 : else
397 0 : ABI_ERROR("Wrong time2")
398 : end if
399 : else
400 0 : if (time2==1) then
401 : !$OMP PARALLEL DO PRIVATE(padat)
402 0 : do dat=1,ndat
403 0 : padat = (dat-1)*nr
404 0 : do ir=1,nr
405 0 : u12prod(ir+padat) = GWPC_CONJG(u1_bz(ir+padat)) * u2_bz(ir+padat)
406 : end do
407 : end do
408 0 : else if (time2==2) then
409 : !$OMP PARALLEL DO PRIVATE(padat)
410 0 : do dat=1,ndat
411 0 : padat = (dat-1)*nr
412 0 : do ir=1,nr
413 0 : u12prod(ir+padat) = GWPC_CONJG(u1_bz(ir+padat)) * GWPC_CONJG(u2_bz(ir+padat))
414 : end do
415 : end do
416 : else
417 0 : ABI_ERROR("Wrong time2")
418 : end if
419 : end if
420 :
421 : CASE (2)
422 367100 : if (ndat==1) then
423 367100 : if (time2==1) then
424 805808192 : do ir=1,nr
425 805808192 : u12prod(ir) = u1_bz(ir) * u2_bz(ir)
426 : end do
427 96672 : else if (time2==2) then
428 299523000 : do ir=1,nr
429 299523000 : u12prod(ir) = u1_bz(ir) * GWPC_CONJG(u2_bz(ir))
430 : end do
431 : else
432 0 : ABI_ERROR("Wrong time2")
433 : end if
434 : else
435 0 : if (time2==1) then
436 : !$OMP PARALLEL DO PRIVATE(padat)
437 0 : do dat=1,ndat
438 0 : padat = (dat-1)*nr
439 0 : do ir=1,nr
440 0 : u12prod(ir+padat) = u1_bz(ir+padat) * u2_bz(ir+padat)
441 : end do
442 : end do
443 0 : else if (time2==2) then
444 : !$OMP PARALLEL DO PRIVATE(padat)
445 0 : do dat=1,ndat
446 0 : padat = (dat-1)*nr
447 0 : do ir=1,nr
448 0 : u12prod(ir+padat) = u1_bz(ir+padat) * GWPC_CONJG(u2_bz(ir+padat))
449 : end do
450 : end do
451 : else
452 0 : ABI_ERROR("Wrong time2")
453 : end if
454 : end if
455 : CASE DEFAULT
456 5643847 : ABI_ERROR("Wrong time1")
457 : END SELECT
458 :
459 5643847 : ABI_FREE(u1_bz)
460 5643847 : ABI_FREE(u2_bz)
461 :
462 5643847 : end subroutine usur_kkp_bz
463 : !!***
464 :
465 : !----------------------------------------------------------------------
466 :
467 : !!****f* m_oscillators/gw_box2gsph
468 : !! NAME
469 : !! gw_box2gsph
470 : !!
471 : !! FUNCTION
472 : !! Trasnfer data from the FFT box to the G-sphere.
473 : !!
474 : !! INPUTS
475 : !! nr=number of FFT grid points
476 : !! ndat=Number of wavefunctions to transform.
477 : !! npw=number of plane waves in the sphere
478 : !! igfftg0(npw)=index of G-G_o in the FFT array for each G in the sphere.
479 : !! iarrbox(nr*ndat)=Input array on the FFT mesh
480 : !!
481 : !! OUTPUT
482 : !! oarrsph(npw*ndat)=output array on the sphere.
483 : !!
484 : !! SOURCE
485 :
486 5643847 : subroutine gw_box2gsph(nr, ndat, npw, igfftg0, iarrbox, oarrsph)
487 :
488 : !Arguments ------------------------------------
489 : !scalars
490 : integer,intent(in) :: nr,ndat,npw
491 : !arrays
492 : integer,intent(in) :: igfftg0(npw)
493 : complex(gwp),intent(in) :: iarrbox(nr*ndat)
494 : complex(gwp),intent(out) :: oarrsph(npw*ndat)
495 :
496 : !Local variables-------------------------------
497 : !scalars
498 : integer :: ig,igfft,dat,pgsp,pfft
499 : ! *************************************************************************
500 :
501 5643847 : if (ndat==1) then
502 327180366 : do ig=1,npw
503 321536519 : igfft=igfftg0(ig)
504 327180366 : if (igfft/=0) then
505 : ! G-G0 belongs to the FFT mesh.
506 321536519 : oarrsph(ig) = iarrbox(igfft)
507 : else
508 : ! Set this component to zero.
509 0 : oarrsph(ig) = czero_gw
510 : end if
511 : end do
512 : else
513 : !$OMP PARALLEL DO PRIVATE(pgsp,pfft,igfft)
514 0 : do dat=1,ndat
515 0 : pgsp = (dat-1)*npw
516 0 : pfft = (dat-1)*nr
517 0 : do ig=1,npw
518 0 : igfft=igfftg0(ig)
519 0 : if (igfft/=0) then
520 : ! G-G0 belongs to the FFT mesh.
521 0 : oarrsph(ig+pgsp) = iarrbox(igfft+pfft)
522 : else
523 : ! Set this component to zero.
524 0 : oarrsph(ig+pgsp) = czero_gw
525 : end if
526 : end do
527 : end do
528 : end if
529 :
530 5643847 : end subroutine gw_box2gsph
531 : !!***
532 :
533 : !----------------------------------------------------------------------
534 :
535 : !!****f* m_oscillators/calc_wfwfg
536 : !! NAME
537 : !! calc_wfwfg
538 : !!
539 : !! FUNCTION
540 : !! Calculate the Fourier transform of the product u_{bk}^*(r) u_{b"k}(r)
541 : !! Return values on the FFT box.
542 : !!
543 : !! INPUTS
544 : !! nspinor=number of spinorial components.
545 : !! spinrot(4)=components of the spinor rotation matrix
546 : !!
547 : !! OUTPUT
548 : !!
549 : !! SOURCE
550 :
551 2842 : subroutine calc_wfwfg(ktabr_k, ktabi_k, spinrot, nr, nspinor, ngfft_gw, wfr_jb, wfr_kb, wfg2_jk)
552 :
553 : !Arguments ------------------------------------
554 : !scalars
555 : integer,intent(in) :: ktabi_k,nr,nspinor
556 : !arrays
557 : integer,intent(in) :: ktabr_k(nr),ngfft_gw(18)
558 : real(dp),intent(in) :: spinrot(4)
559 : complex(gwp),intent(in) :: wfr_jb(nr*nspinor),wfr_kb(nr*nspinor)
560 : complex(gwp),intent(out) :: wfg2_jk(nr*nspinor)
561 :
562 : !Local variables-------------------------------
563 : integer,parameter :: ndat1 = 1, fftcache0 = 0, gpu_option_0 = 0
564 : type(fftbox_plan3_t) :: plan
565 : !arrays
566 2842 : complex(gwp),allocatable :: wfr2_dpcplx(:),ujb_bz(:),ukb_bz(:)
567 : ! *************************************************************************
568 :
569 : ! There is no need to take into account phases arising from non-symmorphic
570 : ! operations since the wavefunctions are evaluated at the same k-point.
571 8526 : ABI_MALLOC(wfr2_dpcplx, (nr * nspinor * ndat1))
572 :
573 2842 : if (nspinor == 1) then
574 5300 : select case (ktabi_k)
575 : case (1)
576 15955530 : wfr2_dpcplx = GWPC_CONJG(wfr_jb(ktabr_k)) * wfr_kb(ktabr_k)
577 : case (2)
578 : ! Conjugate the product if time-reversal is used to reconstruct this k-point
579 540320 : wfr2_dpcplx = wfr_jb(ktabr_k) * GWPC_CONJG(wfr_kb(ktabr_k))
580 : case default
581 2730 : ABI_ERROR(sjoin("Wrong ktabi_k:", itoa(ktabi_k)))
582 : end select
583 :
584 112 : else if (nspinor == 2) then
585 224 : ABI_MALLOC(ujb_bz, (nr * nspinor * ndat1))
586 224 : ABI_MALLOC(ukb_bz, (nr * nspinor * ndat1))
587 : ! Use wfr2_dpcplx as workspace array
588 112 : call rotate_spinor(ktabi_k, ktabr_k, cone, spinrot, nr, nspinor, ndat1, wfr_jb, wfr2_dpcplx, ujb_bz)
589 112 : call rotate_spinor(ktabi_k, ktabr_k, cone, spinrot, nr, nspinor, ndat1, wfr_kb, wfr2_dpcplx, ukb_bz)
590 1792224 : wfr2_dpcplx = GWPC_CONJG(ujb_bz) * ukb_bz
591 112 : ABI_FREE(ujb_bz)
592 112 : ABI_FREE(ukb_bz)
593 :
594 : else
595 0 : ABI_ERROR(sjoin("Wrong nspinor:", itoa(nspinor)))
596 : end if
597 :
598 : ! Transform to Fourier space (result in wfg2_jk)
599 2842 : call plan%init(nspinor, ngfft_gw(1:3), ngfft_gw(1:3), ngfft_gw(7), fftcache0, gpu_option_0)
600 2842 : call plan%execute(wfr2_dpcplx, wfg2_jk, -1, nspinor)
601 2842 : call plan%free()
602 2842 : ABI_FREE(wfr2_dpcplx)
603 :
604 19894 : end subroutine calc_wfwfg
605 : !!***
606 :
607 : !----------------------------------------------------------------------
608 :
609 : !!****f* m_oscillators/sym_rhotwgq0
610 : !! NAME
611 : !! sym_rhotwgq0
612 : !!
613 : !! FUNCTION
614 : !! Symmetrization of the oscillator matrix elements <k-q,b1|exp(-i(q+G).r)|k,b2> in the special case of q=0.
615 : !! The matrix elements in the full BZ is obtained from the matrix elements in the IBZ by
616 : !! rotating the wavefunctions and taking into account time reversal symmetry.
617 : !! strictly speaking the symmetrization can be performed only for non-degenerate states.
618 : !!
619 : !! INPUTS
620 : !! Gsph<gsphere_t>=Info on the G-sphere used to describe wavefunctions and W (the largest one is actually stored).
621 : !! npw=Number of G-vectors
622 : !! dim_rtwg=Number of spin-spin combinations, 1 for collinear spin, 4 is nspinor==2 (TODO NOT CODED)
623 : !! itim_k=2 if time reversal is used to reconstruct the k in the BZ, 1 otherwise.
624 : !! isym_k=The index of the symmetry symrec rotains k_IBZ onto k_BZ.
625 : !! rhxtwg_in(dim_rtwg*npw)=The input matrix elements in the IBZ.
626 : !!
627 : !! OUTPUT
628 : !! rhxtwg_sym(dim_rtwg*npw)=The symmetrized matrix elements in the BZ.
629 : !!
630 : !! NOTES
631 : !! Let M_{G}(k,q) =<k-q,b1|exp(-i(q+G).r)|k,b2>
632 : !! At q ==0, supposing non-degenerate bands, one obtains:
633 : !!
634 : !! 1) M_{ SG}( Sk) = e^{-iSG.t} M_{G} (k)
635 : !! 2) M_{-SG}(-Sk) = e^{+iSG.t} M_{G}^* (k)
636 : !!
637 : !! SOURCE
638 :
639 2084436 : function sym_rhotwgq0(itim_k, isym_k, dim_rtwg, npw, rhxtwg_in, Gsph) result(rhxtwg_sym)
640 :
641 : !Arguments ------------------------------------
642 : !scalars
643 : integer,intent(in) :: npw,dim_rtwg,itim_k,isym_k
644 : type(gsphere_t),intent(in) :: Gsph
645 : !arrays
646 : complex(gwp),intent(in) :: rhxtwg_in(dim_rtwg*npw)
647 : complex(gwp) :: rhxtwg_sym(dim_rtwg*npw)
648 :
649 : !Local variables ------------------------------
650 : !scalars
651 : integer :: ig
652 : !************************************************************************
653 :
654 2084436 : ABI_CHECK(dim_rtwg == 1, "dim_rtwg/=1 not coded")
655 :
656 3537388 : SELECT CASE (isym_k)
657 : CASE (1)
658 : ! Fractional translation associated to E is assumed to be (zero,zero,zero).
659 2084436 : SELECT CASE (itim_k)
660 : CASE (1)
661 : ! Identity, no time-reversal. No symmetrization is needed.
662 62072096 : rhxtwg_sym(:) = rhxtwg_in(:)
663 : CASE (2)
664 : ! Identity + Time-reversal.
665 735072 : do ig=1,npw
666 735072 : rhxtwg_sym( Gsph%rottb(ig,itim_k,isym_k) ) = GWPC_CONJG(rhxtwg_in(ig))
667 : end do
668 : CASE DEFAULT
669 1452952 : ABI_ERROR(sjoin("Wrong value of itim_k:", itoa(itim_k)))
670 : END SELECT
671 :
672 : CASE DEFAULT
673 : ! Rotate wavefunctions.
674 2084436 : SELECT CASE (itim_k)
675 : CASE (1)
676 : ! no time-reversal, only rotation.
677 28247440 : do ig=1,npw
678 28247440 : rhxtwg_sym( Gsph%rottb(ig,itim_k,isym_k) ) = rhxtwg_in(ig) * Gsph%phmSGt(ig,isym_k)
679 : end do
680 : CASE (2)
681 : ! time-reversal + spatial rotation.
682 4905888 : do ig=1,npw
683 4905888 : rhxtwg_sym( Gsph%rottb(ig,itim_k,isym_k) ) = GWPC_CONJG( rhxtwg_in(ig) * Gsph%phmSGt(ig,isym_k) )
684 : end do
685 : CASE DEFAULT
686 631484 : ABI_ERROR(sjoin("Wrong value of itim_k:", itoa(itim_k)))
687 : END SELECT
688 : END SELECT
689 :
690 2084436 : end function sym_rhotwgq0
691 : !!***
692 :
693 : !----------------------------------------------------------------------
694 :
695 : !!****f* m_oscillators/rotate_spinor
696 : !! NAME
697 : !! rotate_spinor
698 : !!
699 : !! FUNCTION
700 : !! Return a spinor in the full BZ from its symmetrical image in the IBZ.
701 : !!
702 : !! INPUTS
703 : !!
704 : !! OUTPUT
705 : !!
706 : !! SOURCE
707 :
708 102400 : subroutine rotate_spinor(itim_kbz, ktabr_kbz, ktabp_kbz, spinrot, nr, nspinor, ndat, ug_ibz, cwork, oug_bz)
709 :
710 : !Arguments ------------------------------------
711 : !scalars
712 : integer,intent(in) :: itim_kbz, nr, nspinor, ndat
713 : complex(dp),intent(in) :: ktabp_kbz
714 : !arrays
715 : integer,intent(in) :: ktabr_kbz(nr)
716 : real(dp),intent(in) :: spinrot(4)
717 : complex(gwp),intent(in) :: ug_ibz(nr*nspinor*ndat)
718 : complex(gwp),intent(out) :: cwork(nr*nspinor*ndat), oug_bz(nr*nspinor*ndat)
719 :
720 : !Local variables ------------------------------
721 : !scalars
722 : integer :: ir,ir1,spad0,ispinor
723 : complex(gwp) :: u1a,u1b
724 : !arrays
725 : complex(dp) :: spinrot_cmat1(2,2)
726 : !************************************************************************
727 :
728 102400 : ABI_CHECK(ndat == 1, "ndat > 1 not coded")
729 102400 : ABI_CHECK(nspinor == 2, "nspinor should be 1")
730 :
731 : ! Step 1: Real-space rotation per spinor component (apply ktabr + phase).
732 : ! This is done BEFORE time-reversal, consistent with cgtk_rotate in m_cgtk.F90:
733 : ! u_{Sk}(r) = e^{-2i\pi k_ibz.(R^{-1}\tau)} u_{k_ibz}(R^{-1}(r-\tau))
734 307200 : do ispinor=1,nspinor
735 204800 : spad0 = (ispinor-1) * nr
736 1638707200 : do ir=1,nr
737 1638400000 : ir1 = ktabr_kbz(ir)
738 1638604800 : oug_bz(ir+spad0) = ug_ibz(ir1+spad0) * ktabp_kbz
739 : end do
740 : end do
741 :
742 : ! Step 2: Apply time-reversal AFTER spatial rotation (consistent with cgtk_rotate).
743 : ! \psi_{-k}^1 = (\psi_k^2)^*
744 : ! \psi_{-k}^2 = -(\psi_k^1)^*
745 102400 : if (itim_kbz == 2) then
746 0 : cwork(1:nr) = GWPC_CONJG(oug_bz(nr+1:2*nr))
747 0 : cwork(nr+1:2*nr) = -GWPC_CONJG(oug_bz(1:nr))
748 0 : oug_bz(:) = cwork(:)
749 : end if
750 :
751 : ! Step 3: Rotation in spinor space using the INVERSE spinor rotation matrix (S^\dagger)
752 : ! since the mapping from IBZ to BZ uses the inverse symmetry operation in real space.
753 102400 : spinrot_cmat1(1,1) = spinrot(1) - j_dpc*spinrot(4)
754 102400 : spinrot_cmat1(1,2) =-spinrot(3) - j_dpc*spinrot(2)
755 102400 : spinrot_cmat1(2,1) = spinrot(3) - j_dpc*spinrot(2)
756 102400 : spinrot_cmat1(2,2) = spinrot(1) + j_dpc*spinrot(4)
757 1638502400 : cwork = oug_bz
758 819302400 : do ir=1,nr
759 819200000 : u1a = cwork(ir); u1b = cwork(ir+nr)
760 819200000 : oug_bz(ir) = spinrot_cmat1(1, 1) * u1a + spinrot_cmat1(1, 2) * u1b
761 819302400 : oug_bz(ir+nr) = spinrot_cmat1(2, 1) * u1a + spinrot_cmat1(2, 2) * u1b
762 : end do
763 :
764 102400 : end subroutine rotate_spinor
765 : !!***
766 :
767 : end module m_oscillators
768 : !!***
|
