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