Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_trafo
8 : use m_judft
9 : use m_glob_tofrom_loc
10 : use m_constants
11 : CONTAINS
12 :
13 0 : SUBROUTINE waveftrafo_symm(cmt_out, z_out, cmt, l_real, z_r, z_c, bandi, ndb, &
14 : nk, iop, atoms, mpdata, hybinp, hybdat, kpts, &
15 : sym, jsp, lapw)
16 :
17 : USE m_constants
18 : USE m_wrapper
19 : USE m_types_mpdata
20 : USE m_types_hybinp
21 : USE m_types_hybdat
22 : USE m_types_sym
23 : USE m_types_kpts
24 : USE m_types_atoms
25 : USE m_types_lapw
26 : USE m_juDFT
27 : IMPLICIT NONE
28 :
29 : TYPE(t_mpdata), INTENT(IN) :: mpdata
30 : TYPE(t_hybinp), INTENT(IN) :: hybinp
31 : TYPE(t_hybdat), INTENT(IN) :: hybdat
32 : TYPE(t_sym), INTENT(IN) :: sym
33 : TYPE(t_kpts), INTENT(IN) :: kpts
34 : TYPE(t_atoms), INTENT(IN) :: atoms
35 : TYPE(t_lapw), INTENT(IN) :: lapw
36 :
37 : ! - scalars -
38 : INTEGER, INTENT(IN) :: nk, jsp, ndb
39 : INTEGER, INTENT(IN) :: bandi, iop
40 :
41 : ! - arrays -
42 : COMPLEX, INTENT(IN) :: cmt(:, :, :)
43 : LOGICAL, INTENT(IN) :: l_real
44 : REAL, INTENT(IN) :: z_r(:, :)
45 : COMPLEX, INTENT(IN) :: z_c(:, :)
46 : COMPLEX, INTENT(INOUT) :: cmt_out(hybdat%maxlmindx, atoms%nat, ndb)
47 : COMPLEX, INTENT(INOUT) :: z_out(lapw%nv(jsp), ndb)
48 :
49 : ! - local -
50 :
51 : ! - scalars -
52 : INTEGER :: iatom, iatom1, iiatom, itype, igpt, igpt1, ieq, iiop
53 : INTEGER :: i, l, n, nn, lm0, lm1, lm2
54 : COMPLEX :: cdum
55 :
56 : ! - arrays -
57 : REAL :: rrot(3, 3), invrrot(3, 3)
58 : INTEGER :: g(3), g1(3)
59 : REAL :: tau1(3), rkpt(3), rkpthlp(3), trans(3)
60 0 : COMPLEX :: cmthlp(2*atoms%lmaxd + 1)
61 : LOGICAL :: trs
62 :
63 0 : if (l_real) THEN
64 0 : rrot = transpose(1.0*sym%mrot(:, :, sym%invtab(iop)))
65 0 : invrrot = transpose(1.0*sym%mrot(:, :, iop))
66 0 : trans = sym%tau(:, iop)
67 : else
68 0 : IF (iop <= sym%nop) THEN
69 0 : trs = .false.
70 0 : rrot = transpose(1.0*sym%mrot(:, :, sym%invtab(iop)))
71 0 : invrrot = transpose(1.0*sym%mrot(:, :, iop))
72 0 : trans = sym%tau(:, iop)
73 : ELSE
74 0 : trs = .true.
75 0 : iiop = iop - sym%nop
76 0 : rrot = -transpose(1.0*sym%mrot(:, :, sym%invtab(iiop)))
77 0 : invrrot = -transpose(1.0*sym%mrot(:, :, iiop))
78 0 : trans = sym%tau(:, iiop)
79 : END IF
80 : end if
81 :
82 0 : rkpt = matmul(rrot, kpts%bkf(:, nk))
83 0 : rkpthlp = rkpt
84 0 : rkpt = kpts%to_first_bz(rkpt)
85 0 : g1 = nint(rkpt - rkpthlp)
86 :
87 : ! MT coefficients
88 0 : cmt_out = 0
89 0 : iatom = 0
90 :
91 0 : DO itype = 1, atoms%ntype
92 0 : iiatom = atoms%firstAtom(itype) - 1
93 0 : DO ieq = 1, atoms%neq(itype)
94 0 : iatom = iatom + 1
95 :
96 0 : iatom1 = hybinp%map(iatom, iop)
97 0 : tau1 = hybinp%tvec(:, iatom, iop)
98 :
99 0 : cdum = exp(-ImagUnit*tpi_const*dot_product(rkpt, tau1))
100 :
101 0 : lm0 = 0
102 0 : DO l = 0, atoms%lmax(itype)
103 0 : nn = mpdata%num_radfun_per_l(l, itype)
104 0 : DO n = 1, nn
105 0 : lm1 = lm0 + n
106 0 : lm2 = lm0 + n + 2*l*nn
107 0 : DO i = 1, ndb
108 0 : if (l_real) THEN
109 : cmt_out(lm1:lm2:nn, iatom1, i) = cdum* &
110 0 : matmul(cmt(bandi + i - 1, lm1:lm2:nn, iatom), &
111 0 : sym%d_wgn(-l:l, -l:l, l, iop))
112 : else
113 0 : IF (trs) THEN
114 0 : cmthlp(:2*l + 1) = CONJG(cmt(bandi + i - 1, lm1:lm2:nn, iatom))
115 : ELSE
116 0 : cmthlp(:2*l + 1) = cmt(bandi + i - 1, lm1:lm2:nn, iatom)
117 : END IF
118 0 : cmt_out(lm1:lm2:nn, iatom1, i) = cdum*matmul(cmthlp(:2*l + 1), sym%d_wgn(-l:l, -l:l, l, iop))
119 : end if
120 : END DO
121 : END DO
122 0 : lm0 = lm2
123 : END DO
124 : END DO
125 : END DO
126 :
127 : ! PW coefficients
128 0 : z_out = 0
129 :
130 0 : DO igpt = 1, lapw%nv(jsp)
131 0 : g = INT(matmul(invrrot, lapw%gvec(:, igpt, jsp) + g1))
132 : !determine number of g
133 0 : igpt1 = 0
134 0 : DO i = 1, lapw%nv(jsp)
135 0 : IF (maxval(abs(g - lapw%gvec(:, i, jsp))) <= 1E-06) THEN
136 : igpt1 = i
137 : EXIT
138 : END IF
139 : END DO
140 0 : IF (igpt1 == 0) THEN
141 0 : call judft_error('wavetrafo_symm: rotated G vector not found')
142 : END IF
143 0 : cdum = exp(-ImagUnit*tpi_const*dot_product(rkpt + lapw%gvec(:, igpt, jsp), trans(:)))
144 0 : if (l_real) THEN
145 0 : z_out(igpt, 1:ndb) = cdum*z_r(igpt1, bandi:bandi + ndb - 1)
146 : else
147 0 : IF (trs) THEN
148 0 : z_out(igpt, 1:ndb) = cdum*CONJG(z_c(igpt1, bandi:bandi + ndb - 1))
149 : ELSE
150 0 : z_out(igpt, 1:ndb) = cdum*z_c(igpt1, bandi:bandi + ndb - 1)
151 : END IF
152 : end if
153 : END DO
154 :
155 0 : END SUBROUTINE waveftrafo_symm
156 :
157 84 : SUBROUTINE waveftrafo_gen_cmt(cmt, c_phase, l_real, nk, iop, atoms, &
158 28 : mpdata, hybinp, kpts, sym, nbands, cmt_out)
159 :
160 : use m_juDFT
161 : USE m_constants
162 : USE m_wrapper
163 : USE m_types_mpdata
164 : USE m_types_hybinp
165 : USE m_types_sym
166 : USE m_types_kpts
167 : USE m_types_atoms
168 : IMPLICIT NONE
169 :
170 : TYPE(t_mpdata), INTENT(IN) :: mpdata
171 : TYPE(t_hybinp), INTENT(IN) :: hybinp
172 : TYPE(t_sym), INTENT(IN) :: sym
173 : TYPE(t_kpts), INTENT(IN) :: kpts
174 : TYPE(t_atoms), INTENT(IN) :: atoms
175 : ! - scalars -
176 : INTEGER, INTENT(IN) :: nk, nbands
177 : INTEGER, INTENT(IN) :: iop
178 : LOGICAL, INTENT(in) :: l_real
179 : ! - arrays -
180 : COMPLEX, INTENT(IN) :: cmt(:, :, :), c_phase(nbands)
181 :
182 : COMPLEX, INTENT(INOUT) :: cmt_out(:, :, :)
183 : ! - local -
184 :
185 : ! - scalars -
186 : INTEGER :: itype, iatom, iatom1, iiatom, ieq, iiop
187 : INTEGER :: i, l, n, nn, lm0, lm1, lm2
188 : COMPLEX :: cdum
189 : LOGICAL :: trs
190 :
191 : ! - arrays -
192 : INTEGER :: rrot(3, 3), invrrot(3, 3)
193 : INTEGER :: g1(3)
194 : REAL :: tau1(3), rkpt(3), rkpthlp(3), trans(3)
195 28 : COMPLEX :: cmthlp(2*atoms%lmaxd + 1)
196 :
197 28 : call timestart("gen_cmt")
198 28 : if (l_real) THEN
199 260 : rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
200 : invrrot = transpose(sym%mrot(:, :, iop))
201 20 : trans = sym%tau(:, iop)
202 : else
203 8 : IF (iop <= sym%nop) THEN
204 8 : trs = .false.
205 104 : rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
206 : invrrot = transpose(sym%mrot(:, :, iop))
207 8 : trans = sym%tau(:, iop)
208 : ELSE
209 : ! in the case of SOC (l_soc=.true.)
210 : ! time reversal symmetry is not valid anymore;
211 : ! nsym should thus equal nop
212 0 : trs = .true.
213 0 : iiop = iop - sym%nop
214 0 : rrot = -transpose(sym%mrot(:, :, sym%invtab(iiop)))
215 : invrrot = -transpose(sym%mrot(:, :, iiop))
216 0 : trans = sym%tau(:, iiop)
217 : END IF
218 : end if
219 :
220 700 : rkpt = matmul(rrot, kpts%bkf(:, nk))
221 : rkpthlp = rkpt
222 112 : rkpt = kpts%to_first_bz(rkpt)
223 28 : g1 = nint(rkpt - rkpthlp)
224 :
225 : ! MT coefficients
226 93504 : cmt_out = 0
227 28 : iatom = 0
228 :
229 72 : DO itype = 1, atoms%ntype
230 44 : iiatom = atoms%firstAtom(itype) - 1
231 116 : DO ieq = 1, atoms%neq(itype)
232 44 : iatom = iatom + 1
233 :
234 44 : iatom1 = hybinp%map(iatom, iop)
235 176 : tau1 = hybinp%tvec(:, iatom, iop)
236 :
237 176 : cdum = exp(-ImagUnit*tpi_const*dot_product(rkpt, tau1))
238 :
239 44 : lm0 = 0
240 516 : DO l = 0, atoms%lmax(itype)
241 428 : nn = mpdata%num_radfun_per_l(l, itype)
242 1340 : DO n = 1, nn
243 912 : lm1 = lm0 + n
244 912 : lm2 = lm0 + n + 2*l*nn
245 :
246 10460 : DO i = 1, nbands
247 10032 : if (l_real) THEN
248 4864 : cmt_out(i, lm1:lm2:nn, iatom1) = cdum*matmul(cmt(i, lm1:lm2:nn, iatom), &
249 787392 : hybinp%d_wgn2(-l:l, -l:l, l, iop))
250 : else
251 4256 : IF (trs) THEN
252 0 : cmthlp(:2*l + 1) = conjg(cmt(i, lm1:lm2:nn, iatom))
253 : ELSE
254 41664 : cmthlp(:2*l + 1) = cmt(i, lm1:lm2:nn, iatom)
255 : END IF
256 556192 : cmt_out(i, lm1:lm2:nn, iatom1) = cdum*matmul(cmthlp(:2*l + 1), hybinp%d_wgn2(-l:l, -l:l, l, iop))
257 : end if
258 :
259 : END DO
260 : END DO
261 472 : lm0 = lm2
262 : END DO
263 : END DO
264 : END DO
265 :
266 : ! If phase and inversion-sym. is true,
267 : ! define the phase such that z_out is real.
268 28 : if (l_real) then
269 180 : DO i = 1, nbands
270 47828 : cmt_out(i, :, :) = cmt_out(i, :, :)/c_phase(i)
271 : END DO
272 : end if
273 28 : call timestop("gen_cmt")
274 28 : END SUBROUTINE waveftrafo_gen_cmt
275 :
276 0 : SUBROUTINE waveftrafo_genwavf( &
277 0 : cmt, z_in, nk, iop, atoms, &
278 : mpdata, hybinp, kpts, sym, jsp, input, nbands, &
279 0 : lapw_nk, lapw_rkpt, cmt_out, z_out)
280 :
281 : use m_juDFT
282 : USE m_constants
283 : USE m_wrapper
284 : USE m_types_mat
285 : USE m_types_input
286 : USE m_types_mpdata
287 : USE m_types_hybinp
288 : USE m_types_sym
289 : USE m_types_kpts
290 : USE m_types_atoms
291 : USE m_types_lapw
292 : IMPLICIT NONE
293 :
294 : type(t_mat), intent(in) :: z_in
295 : TYPE(t_input), INTENT(IN) :: input
296 : TYPE(t_mpdata), INTENT(IN) :: mpdata
297 : TYPE(t_hybinp), INTENT(IN) :: hybinp
298 : TYPE(t_sym), INTENT(IN) :: sym
299 : TYPE(t_kpts), INTENT(IN) :: kpts
300 : TYPE(t_atoms), INTENT(IN) :: atoms
301 : TYPE(t_lapw), INTENT(IN) :: lapw_nk, lapw_rkpt
302 : type(t_mat), intent(inout) :: z_out
303 : ! - scalars -
304 : INTEGER, INTENT(IN) :: nk, jsp, nbands
305 : INTEGER, INTENT(IN) :: iop
306 : ! - arrays -
307 : COMPLEX, INTENT(IN) :: cmt(:, :, :)
308 :
309 : COMPLEX, INTENT(INOUT) :: cmt_out(:, :, :)
310 : ! - local -
311 :
312 : ! - scalars -
313 : INTEGER :: itype, iatom, iatom1, iiatom, igpt, igpt1, ieq, iiop
314 : INTEGER :: i, l, n, nn, lm0, lm1, lm2
315 : COMPLEX :: cdum
316 : LOGICAL :: trs
317 :
318 : ! - arrays -
319 : INTEGER :: rrot(3, 3), invrrot(3, 3)
320 : INTEGER :: g(3), g1(3)
321 : REAL :: tau1(3), rkpt(3), rkpthlp(3), trans(3)
322 0 : COMPLEX :: zhlp(z_in%matsize1, input%neig)
323 0 : COMPLEX :: cmthlp(2*atoms%lmaxd + 1)
324 :
325 0 : call timestart("genwavf")
326 0 : if (z_in%l_real) THEN
327 0 : rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
328 0 : invrrot = transpose(sym%mrot(:, :, iop))
329 0 : trans = sym%tau(:, iop)
330 : else
331 0 : IF (iop <= sym%nop) THEN
332 0 : trs = .false.
333 0 : rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
334 0 : invrrot = transpose(sym%mrot(:, :, iop))
335 0 : trans = sym%tau(:, iop)
336 : ELSE
337 : ! in the case of SOC (l_soc=.true.)
338 : ! time reversal symmetry is not valid anymore;
339 : ! nsym should thus equal nop
340 0 : trs = .true.
341 0 : iiop = iop - sym%nop
342 0 : rrot = -transpose(sym%mrot(:, :, sym%invtab(iiop)))
343 0 : invrrot = -transpose(sym%mrot(:, :, iiop))
344 0 : trans = sym%tau(:, iiop)
345 : END IF
346 : end if
347 :
348 0 : rkpt = matmul(rrot, kpts%bkf(:, nk))
349 0 : rkpthlp = rkpt
350 0 : rkpt = kpts%to_first_bz(rkpt)
351 0 : g1 = nint(rkpt - rkpthlp)
352 :
353 : ! MT coefficients
354 0 : cmt_out = 0
355 0 : iatom = 0
356 :
357 0 : DO itype = 1, atoms%ntype
358 0 : iiatom = atoms%firstAtom(itype) - 1
359 0 : DO ieq = 1, atoms%neq(itype)
360 0 : iatom = iatom + 1
361 :
362 0 : iatom1 = hybinp%map(iatom, iop)
363 0 : tau1 = hybinp%tvec(:, iatom, iop)
364 :
365 0 : cdum = exp(-ImagUnit*tpi_const*dot_product(rkpt, tau1))
366 :
367 0 : lm0 = 0
368 0 : DO l = 0, atoms%lmax(itype)
369 0 : nn = mpdata%num_radfun_per_l(l, itype)
370 0 : DO n = 1, nn
371 0 : lm1 = lm0 + n
372 0 : lm2 = lm0 + n + 2*l*nn
373 :
374 0 : DO i = 1, nbands
375 0 : if (z_in%l_real) THEN
376 0 : cmt_out(i, lm1:lm2:nn, iatom1) = cdum*matmul(cmt(i, lm1:lm2:nn, iatom), &
377 0 : hybinp%d_wgn2(-l:l, -l:l, l, iop))
378 : else
379 0 : IF (trs) THEN
380 0 : cmthlp(:2*l + 1) = conjg(cmt(i, lm1:lm2:nn, iatom))
381 : ELSE
382 0 : cmthlp(:2*l + 1) = cmt(i, lm1:lm2:nn, iatom)
383 : END IF
384 0 : cmt_out(i, lm1:lm2:nn, iatom1) = cdum*matmul(cmthlp(:2*l + 1), hybinp%d_wgn2(-l:l, -l:l, l, iop))
385 : end if
386 :
387 : END DO
388 : END DO
389 0 : lm0 = lm2
390 : END DO
391 : END DO
392 : END DO
393 :
394 : ! PW coefficients
395 :
396 0 : zhlp = 0
397 0 : DO igpt = 1, lapw_rkpt%nv(jsp)
398 0 : g = matmul(invrrot, lapw_rkpt%gvec(:, igpt, jsp) + g1)
399 : !determine number of g
400 0 : igpt1 = 0
401 0 : DO i = 1, lapw_nk%nv(jsp)
402 0 : IF (all(abs(g - lapw_nk%gvec(:, i, jsp)) <= 1E-06)) THEN
403 : igpt1 = i
404 : EXIT
405 : END IF
406 : END DO
407 0 : IF (igpt1 == 0) CYCLE
408 0 : cdum = exp(-ImagUnit*tpi_const*dot_product(rkpt + lapw_rkpt%gvec(:, igpt, jsp), trans))
409 0 : if (z_in%l_real) THEN
410 0 : zhlp(igpt, :nbands) = cdum*z_in%data_r(igpt1, :nbands)
411 : else
412 0 : IF (trs) THEN
413 0 : zhlp(igpt, :nbands) = cdum*conjg(z_in%data_c(igpt1, :nbands))
414 : ELSE
415 0 : zhlp(igpt, :nbands) = cdum*z_in%data_c(igpt1, :nbands)
416 : END IF
417 : end if
418 : END DO
419 :
420 : ! If phase and inversion-sym. is true,
421 : ! define the phase such that z_out is real.
422 :
423 0 : DO i = 1, nbands
424 0 : if (z_in%l_real) THEN
425 0 : cdum = commonphase(zhlp(:, i), z_in%matsize1)
426 :
427 0 : IF (any(abs(aimag(zhlp(:, i)/cdum)) > 1e-8)) THEN
428 0 : WRITE (*, *) maxval(abs(aimag(zhlp(:, i)/cdum)))
429 0 : WRITE (*, *) zhlp
430 0 : call judft_error('waveftrafo1: Residual imaginary part.')
431 : END IF
432 0 : z_out%data_r(:, i) = real(zhlp(:, i)/cdum)
433 0 : cmt_out(i, :, :) = cmt_out(i, :, :)/cdum
434 : else
435 0 : z_out%data_c(:, i) = zhlp(:, i)
436 : end if
437 : END DO
438 0 : call timestop("genwavf")
439 0 : END SUBROUTINE waveftrafo_genwavf
440 :
441 28 : SUBROUTINE waveftrafo_gen_zmat(z_in, nk, iop, &
442 : kpts, sym, jsp, nbands, &
443 28 : lapw_nk, lapw_rkpt, z_out, c_phase)
444 :
445 : use m_juDFT
446 : USE m_constants
447 : USE m_wrapper
448 : USE m_types_mat
449 : USE m_types_sym
450 : USE m_types_kpts
451 : USE m_types_lapw
452 : IMPLICIT NONE
453 :
454 : type(t_mat), intent(in) :: z_in
455 : TYPE(t_sym), INTENT(IN) :: sym
456 : TYPE(t_kpts), INTENT(IN) :: kpts
457 : TYPE(t_lapw), INTENT(IN) :: lapw_nk, lapw_rkpt
458 : type(t_mat), intent(inout) :: z_out
459 : complex, intent(inout), optional :: c_phase(:)
460 : ! - scalars -
461 : INTEGER, INTENT(IN) :: nk, jsp, nbands
462 : INTEGER, INTENT(IN) :: iop
463 :
464 : ! - scalars -
465 : INTEGER :: igpt, igpt1, iiop, i
466 : COMPLEX :: cdum
467 : LOGICAL :: trs
468 :
469 : ! - arrays -
470 : INTEGER :: rrot(3, 3), invrrot(3, 3)
471 : INTEGER :: g(3), g1(3)
472 : REAL :: rkpt(3), rkpthlp(3), trans(3)
473 28 : COMPLEX :: zhlp(z_in%matsize1, nbands)
474 :
475 28 : call timestart("gen_zmat")
476 1440 : if (present(c_phase)) c_phase = 0
477 :
478 28 : if (z_in%l_real) THEN
479 260 : rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
480 260 : invrrot = transpose(sym%mrot(:, :, iop))
481 80 : trans = sym%tau(:, iop)
482 : else
483 8 : IF (iop <= sym%nop) THEN
484 8 : trs = .false.
485 104 : rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
486 104 : invrrot = transpose(sym%mrot(:, :, iop))
487 32 : trans = sym%tau(:, iop)
488 : ELSE
489 : ! in the case of SOC (l_soc=.true.)
490 : ! time reversal symmetry is not valid anymore;
491 : ! nsym should thus equal nop
492 0 : trs = .true.
493 0 : iiop = iop - sym%nop
494 0 : rrot = -transpose(sym%mrot(:, :, sym%invtab(iiop)))
495 0 : invrrot = -transpose(sym%mrot(:, :, iiop))
496 0 : trans = sym%tau(:, iiop)
497 : END IF
498 : end if
499 :
500 700 : rkpt = matmul(rrot, kpts%bkf(:, nk))
501 28 : rkpthlp = rkpt
502 112 : rkpt = kpts%to_first_bz(rkpt)
503 112 : g1 = nint(rkpt - rkpthlp)
504 :
505 : ! PW coefficients
506 :
507 43612 : zhlp = 0
508 4196 : DO igpt = 1, lapw_rkpt%nv(jsp)
509 66688 : g = matmul(invrrot, lapw_rkpt%gvec(:, igpt, jsp) + g1)
510 : !determine number of g
511 4168 : igpt1 = 0
512 376668 : DO i = 1, lapw_nk%nv(jsp)
513 452200 : IF (all(abs(g - lapw_nk%gvec(:, i, jsp)) <= 1E-06)) THEN
514 : igpt1 = i
515 : EXIT
516 : END IF
517 : END DO
518 4168 : IF (igpt1 == 0) CYCLE
519 16672 : cdum = exp(-ImagUnit*tpi_const*dot_product(rkpt + lapw_rkpt%gvec(:, igpt, jsp), trans))
520 4196 : if (z_in%l_real) THEN
521 25200 : zhlp(igpt, :nbands) = cdum*z_in%data_r(igpt1, :nbands)
522 : else
523 1368 : IF (trs) THEN
524 0 : zhlp(igpt, :nbands) = cdum*conjg(z_in%data_c(igpt1, :nbands))
525 : ELSE
526 20520 : zhlp(igpt, :nbands) = cdum*z_in%data_c(igpt1, :nbands)
527 : END IF
528 : end if
529 : END DO
530 :
531 : ! If phase and inversion-sym. is true,
532 : ! define the phase such that z_out is real.
533 :
534 300 : DO i = 1, nbands
535 300 : if (z_in%l_real) THEN
536 160 : cdum = commonphase(zhlp(:, i), z_in%matsize1)
537 160 : if (present(c_phase)) c_phase(i) = cdum
538 160 : if (abs(cdum) < 1e-30) THEN
539 0 : call juDFT_error("commonphase can't be 0.")
540 : end if
541 23200 : IF (any(abs(aimag(zhlp(:, i)/cdum)) > 1e-8)) THEN
542 0 : WRITE (*, *) maxval(abs(aimag(zhlp(:, i)/cdum)))
543 0 : call judft_error('waveftrafo1: Residual imaginary part.')
544 : END IF
545 23200 : z_out%data_r(:, i) = real(zhlp(:, i)/cdum)
546 : else
547 20384 : z_out%data_c(:, i) = zhlp(:, i)
548 : end if
549 : END DO
550 28 : call timestop("gen_zmat")
551 28 : END SUBROUTINE waveftrafo_gen_zmat
552 :
553 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
554 : ! Symmetrizes MT part of input matrix according to inversion symmetry.
555 : ! This is achieved by a transformation to
556 : ! 1/sqrt(2) * ( exp(ikR) Y_lm(r-R) + (-1)**(l+m) exp(-ikR) Y_l,-m(r+R) )
557 : ! and if R /=0 or m<0
558 : ! i/sqrt(2) * ( exp(ikR) Y_lm(r-R) - (-1)**(l+m) exp(-ikR) Y_l,-m(r+R) ) .
559 : !
560 : ! or
561 : ! i*Y_l,0(r) if R=0,m=0 and l odd
562 : ! These functions have the property f(-r)=f(r)* which makes the output matrix real symmetric.
563 : ! (Array mat is overwritten! )
564 :
565 120 : SUBROUTINE symmetrize_mpimat(fi, fmpi, mpimat, start_dim, end_dim, imode, lreal, nindxm)
566 : USE m_types_fleurinput
567 : USE m_types_mpi
568 : use m_constants
569 :
570 : #ifdef CPP_MPI
571 : USE mpi
572 : #endif
573 :
574 :
575 : IMPLICIT NONE
576 : type(t_fleurinput), intent(in) :: fi
577 : type(t_mpi), intent(in) :: fmpi
578 :
579 : ! - scalars -
580 : INTEGER, INTENT(IN) :: imode, start_dim(2), end_dim(2)
581 : LOGICAL, INTENT(IN) :: lreal
582 :
583 : ! - arrays -
584 : INTEGER, INTENT(IN) :: nindxm(0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype)
585 : COMPLEX, INTENT(INOUT) :: mpimat(:, :)
586 :
587 : ! -local scalars -
588 : INTEGER :: i, j, itype, ieq, ic, ic1, l, m, n, nn, ifac, ishift, start_loc, end_loc
589 : integer :: i_loc, j_loc, i_pe, j_pe, ierr, len_dim(2)
590 : REAL :: rfac
591 :
592 : ! - local arrays -
593 144 : COMPLEX :: mpicarr(maxval(end_dim)), cfac
594 :
595 48 : call timestart("symmetrize_mpimat")
596 144 : len_dim = end_dim - start_dim + 1
597 48 : rfac = sqrt(0.5)
598 48 : cfac = sqrt(0.5)*ImagUnit
599 48 : ic = 0
600 48 : i = 0
601 :
602 120 : DO itype = 1, fi%atoms%ntype
603 864 : nn = sum([((2*l + 1)*nindxm(l, itype), l=0, fi%hybinp%lcutm1(itype))])
604 192 : DO ieq = 1, fi%atoms%neq(itype)
605 72 : ic = ic + 1
606 72 : IF (fi%sym%invsat(ic) == 0) THEN
607 : ! if the structure is inversion-symmetric, but the equivalent atom belongs to a different unit cell
608 : ! invsat(atom) = 0, invsatnr(atom) = 0
609 : ! but we need invsatnr(atom) = natom
610 : ic1 = ic
611 : ELSE
612 0 : ic1 = fi%sym%invsatnr(ic)
613 : END IF
614 : !ic1 = invsatnr(ic)
615 72 : IF (ic1 < ic) THEN
616 0 : i = i + nn
617 0 : CYCLE
618 : END IF
619 : ! IF( ic1 .lt. ic ) cycle
620 504 : DO l = 0, fi%hybinp%lcutm1(itype)
621 360 : ifac = -1
622 2232 : DO m = -l, l
623 1800 : ifac = -ifac
624 1800 : ishift = (ic1 - ic)*nn - 2*m*nindxm(l, itype)
625 14448 : DO n = 1, nindxm(l, itype)
626 12288 : i = i + 1
627 12288 : j = i + ishift
628 12288 : call glob_to_loc(fmpi, i, i_pe, i_loc)
629 12288 : call glob_to_loc(fmpi, j, j_pe, j_loc)
630 14088 : IF (ic1 /= ic .or. m < 0) THEN
631 4800 : IF (iand(imode, 1) /= 0) THEN
632 4800 : call range_from_glob_to_loc(fmpi, start_dim(2), start_loc)
633 4800 : call range_to_glob_to_loc(fmpi, end_dim(2), end_loc)
634 513088 : mpicarr(start_loc:end_loc) = mpimat(i,start_loc:end_loc)
635 513088 : mpimat(i,start_loc:end_loc) = (mpicarr(start_loc:end_loc) + ifac*mpimat(j,start_loc:end_loc))*rfac
636 513088 : mpimat(j,start_loc:end_loc) = (mpicarr(start_loc:end_loc) - ifac*mpimat(j,start_loc:end_loc))*(-cfac)
637 : END IF
638 4800 : IF (iand(imode, 2) /= 0) THEN
639 2400 : if(i_pe == j_pe .and. fmpi%n_rank == i_pe) then
640 332256 : mpicarr(start_dim(1):end_dim(1)) = mpimat(start_dim(1):end_dim(1), i_loc)
641 : mpimat(start_dim(1):end_dim(1),i_loc) &
642 332256 : = (mpimat(start_dim(1):end_dim(1), i_loc) + ifac*mpimat(start_dim(1):end_dim(1), j_loc))*rfac
643 : mpimat(start_dim(1):end_dim(1),j_loc) &
644 332256 : = (mpicarr(start_dim(1):end_dim(1)) - ifac*mpimat(start_dim(1):end_dim(1), j_loc))*cfac
645 : #ifdef CPP_MPI
646 : else
647 1200 : if(fmpi%n_rank == i_pe) then
648 0 : call MPI_Send(mpimat(start_dim(1),i_loc), len_dim(1), MPI_DOUBLE_COMPLEX, j_pe, i, fmpi%sub_comm, ierr)
649 0 : call MPI_Recv(mpicarr(start_dim(1)), len_dim(1), MPI_DOUBLE_COMPLEX, j_pe, j, fmpi%sub_comm, MPI_STATUS_IGNORE, ierr)
650 0 : mpimat(start_dim(1):end_dim(1),i_loc) = (mpimat(start_dim(1):end_dim(1), i_loc) + ifac*mpicarr(start_dim(1):end_dim(1)))*rfac
651 1200 : elseif(fmpi%n_rank == j_pe) then
652 0 : call MPI_Recv(mpicarr(start_dim(1)), len_dim(1), MPI_DOUBLE_COMPLEX, i_pe, i, fmpi%sub_comm, MPI_STATUS_IGNORE, ierr)
653 0 : call MPI_Send(mpimat(start_dim(1),j_loc), len_dim(1), MPI_DOUBLE_COMPLEX, i_pe, j, fmpi%sub_comm, ierr)
654 0 : mpimat(start_dim(1):end_dim(1),j_loc) = (mpicarr(start_dim(1):end_dim(1)) - ifac*mpimat(start_dim(1):end_dim(1), j_loc))*cfac
655 : endif
656 : #endif
657 : endif
658 : END IF
659 7488 : ELSE IF (m == 0 .and. ifac == -1) THEN
660 1056 : IF (iand(imode, 1) /= 0) THEN
661 1056 : call range_from_glob_to_loc(fmpi, start_dim(2), start_loc)
662 1056 : call range_to_glob_to_loc(fmpi, end_dim(2), end_loc)
663 112592 : mpimat(i,start_loc:end_loc) = -ImagUnit*mpimat(i,start_loc:end_loc)
664 : END IF
665 1056 : IF (iand(imode, 2) /= 0 .and. fmpi%n_rank == i_pe) THEN
666 72960 : mpimat(start_dim(1):end_dim(1), i_loc) = ImagUnit*mpimat(start_dim(1):end_dim(1), i_loc)
667 : END IF
668 : END IF
669 : END DO
670 : END DO
671 : END DO
672 : END DO
673 : END DO
674 :
675 48 : IF (lreal) THEN
676 0 : call juDFT_error("this isn't impemented for mpimat")
677 : ! Determine common phase factor and divide by it to make the output matrix real.
678 : ! cfac = commonphase_mtx(mat, dims(1), dims(2))
679 : ! do i = 1, dims(1)
680 : ! do j = 1, dims(2)
681 : ! mat(i, j) = mat(i, j)/cfac
682 : ! if (abs(aimag(mat(i, j))) > 1e-8) then
683 : ! call judft_error('symmetrize: Residual imaginary part. Symmetrization failed.')
684 : ! end if
685 : ! end do
686 : ! end do
687 : END IF
688 :
689 48 : call timestop("symmetrize_mpimat")
690 48 : END SUBROUTINE symmetrize_mpimat
691 :
692 :
693 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
694 : ! Symmetrizes MT part of input matrix according to inversion symmetry.
695 : ! This is achieved by a transformation to
696 : ! 1/sqrt(2) * ( exp(ikR) Y_lm(r-R) + (-1)**(l+m) exp(-ikR) Y_l,-m(r+R) )
697 : ! and if R /=0 or m<0
698 : ! i/sqrt(2) * ( exp(ikR) Y_lm(r-R) - (-1)**(l+m) exp(-ikR) Y_l,-m(r+R) ) .
699 : !
700 : ! or
701 : ! i*Y_l,0(r) if R=0,m=0 and l odd
702 : ! These functions have the property f(-r)=f(r)* which makes the output matrix real symmetric.
703 : ! (Array mat is overwritten! )
704 :
705 22152 : SUBROUTINE symmetrize(mat, dim1, dim2, imode,&
706 7384 : atoms, lcutm, maxlcutm, nindxm, sym)
707 : USE m_types_atoms
708 : USE m_types_sym
709 : use m_constants
710 : IMPLICIT NONE
711 : TYPE(t_atoms), INTENT(IN) :: atoms
712 : TYPE(t_sym), INTENT(IN) :: sym
713 :
714 : ! - scalars -
715 : INTEGER, INTENT(IN) :: imode, dim1, dim2
716 : INTEGER, INTENT(IN) :: maxlcutm
717 :
718 : ! - arrays -
719 : INTEGER, INTENT(IN) :: lcutm(:)
720 : INTEGER, INTENT(IN) :: nindxm(0:maxlcutm, atoms%ntype)
721 : COMPLEX, INTENT(INOUT) :: mat(:, :)
722 :
723 : ! -local scalars -
724 : INTEGER :: i, j, itype, ieq, ic, ic1, l, m, n, nn, ifac, ishift
725 : REAL, parameter :: rfac = sqrt(0.5)
726 :
727 : ! - local arrays -
728 7384 : COMPLEX :: carr(max(dim1, dim2)), cfac = sqrt(0.5)*ImagUnit
729 :
730 : ! call timestart("symmetrize")
731 7384 : ic = 0
732 7384 : i = 0
733 :
734 18030 : DO itype = 1, atoms%ntype
735 127608 : nn = sum([((2*l + 1)*nindxm(l, itype), l=0, lcutm(itype))])
736 28676 : DO ieq = 1, atoms%neq(itype)
737 10646 : ic = ic + 1
738 10646 : IF (sym%invsat(ic) == 0) THEN
739 : ! if the structure is inversion-symmetric, but the equivalent atom belongs to a different unit cell
740 : ! invsat(atom) = 0, invsatnr(atom) = 0
741 : ! but we need invsatnr(atom) = natom
742 : ic1 = ic
743 : ELSE
744 0 : ic1 = sym%invsatnr(ic)
745 : END IF
746 : !ic1 = invsatnr(ic)
747 10646 : IF (ic1 < ic) THEN
748 0 : i = i + nn
749 0 : CYCLE
750 : END IF
751 : ! IF( ic1 .lt. ic ) cycle
752 74450 : DO l = 0, lcutm(itype)
753 53158 : ifac = -1
754 329450 : DO m = -l, l
755 265646 : ifac = -ifac
756 265646 : ishift = (ic1 - ic)*nn - 2*m*nindxm(l, itype)
757 2143986 : DO n = 1, nindxm(l, itype)
758 1825182 : i = i + 1
759 1825182 : j = i + ishift
760 2090828 : IF (ic1 /= ic .or. m < 0) THEN
761 712712 : IF (iand(imode, 1) /= 0) THEN
762 1417580 : carr(:dim2) = mat(i, :dim2)
763 1417580 : mat(i, :dim2) = (carr(:dim2) + ifac*mat(j, :dim2))*rfac
764 1417580 : mat(j, :dim2) = (carr(:dim2) - ifac*mat(j, :dim2))*(-cfac)
765 : END IF
766 712712 : IF (iand(imode, 2) /= 0) THEN
767 8204 : carr(:dim1) = mat(:dim1, i)
768 8204 : mat(:dim1, i) = (carr(:dim1) + ifac*mat(:dim1, j))*rfac
769 8204 : mat(:dim1, j) = (carr(:dim1) - ifac*mat(:dim1, j))*cfac
770 : END IF
771 1112470 : ELSE IF (m == 0 .and. ifac == -1) THEN
772 156952 : IF (iand(imode, 1) /= 0) THEN
773 312300 : mat(i, :dim2) = -ImagUnit*mat(i, :dim2)
774 : END IF
775 156952 : IF (iand(imode, 2) /= 0) THEN
776 1964 : mat(:dim1, i) = ImagUnit*mat(:dim1, i)
777 : END IF
778 : END IF
779 : END DO
780 : END DO
781 : END DO
782 : END DO
783 : END DO
784 : ! call timestop("symmetrize")
785 7384 : END SUBROUTINE symmetrize
786 :
787 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
788 :
789 : ! Undoes symmetrization with routine symmetrize.
790 14652 : SUBROUTINE desymmetrize(mat, dim1, dim2, &
791 7326 : atoms, lcutm, maxlcutm, nindxm, sym)
792 :
793 : USE m_types_sym
794 : USE m_types_atoms
795 : IMPLICIT NONE
796 : TYPE(t_sym), INTENT(IN) :: sym
797 : TYPE(t_atoms), INTENT(IN) :: atoms
798 :
799 : ! - scalars -
800 : INTEGER, INTENT(IN) :: dim1, dim2
801 : INTEGER, INTENT(IN) :: maxlcutm
802 :
803 : ! - arrays -
804 : INTEGER, INTENT(IN) :: lcutm(:)
805 : INTEGER, INTENT(IN) :: nindxm(0:maxlcutm, atoms%ntype)
806 : COMPLEX, INTENT(INOUT) :: mat(dim1, dim2)
807 :
808 : ! - local scalars -
809 : INTEGER :: ifac, i, istart, j, itype, ieq, ic, ic1, l, m, n, nn, ishift
810 : REAL, parameter :: rfac1 = sqrt(0.5)
811 : real :: rfac2
812 : ! - local arrays -
813 7326 : COMPLEX :: carr(max(dim1, dim2))
814 :
815 : ! call timestart("desymmetrize")
816 7326 : ic = 0
817 7326 : istart = 0
818 17888 : DO itype = 1, atoms%ntype
819 126744 : nn = sum([((2*l + 1)*nindxm(l, itype), l=0, lcutm(itype))])
820 28450 : DO ieq = 1, atoms%neq(itype)
821 10562 : ic = ic + 1
822 : ! if the structure is inversion-symmetric, but the equivalent atom belongs to a different unit cell
823 : ! invsat(atom) = 0, invsatnr(atom) =0
824 : ! but we need invsatnr(atom) = natom
825 10562 : ic1 = merge(ic, sym%invsatnr(ic), sym%invsat(ic) == 0)
826 : !ic1 = invsatnr(ic)
827 : !IF( ic1 .lt. ic ) cycle
828 21124 : IF (ic1 < ic) THEN
829 0 : istart = istart + nn
830 : else
831 63372 : DO l = 0, lcutm(itype)
832 52810 : ifac = -1
833 327422 : DO m = -l, l
834 264050 : ifac = -ifac
835 264050 : rfac2 = rfac1*ifac
836 264050 : ishift = (ic1 - ic)*nn - 2*m*nindxm(l, itype)
837 264050 : IF (ic1 /= ic .or. m < 0) THEN
838 105620 : if (ishift <= nindxm(l, itype)) call juDFT_error("if ishift is zero the parallelization is wrong")
839 814308 : DO n = 1, nindxm(l, itype)
840 708688 : i = istart + n
841 708688 : j = i + ishift
842 1417376 : carr(:dim2) = mat(i, :)
843 1417376 : mat(i, :) = (carr(:dim2) + ImagUnit*mat(j, :))*rfac1
844 1522996 : mat(j, :) = (carr(:dim2) - ImagUnit*mat(j, :))*rfac2
845 : enddo
846 158430 : ELSE IF (m == 0 .and. ifac == -1) THEN
847 177172 : DO n = 1, nindxm(l, itype)
848 333220 : mat(istart + n, :) = ImagUnit*mat(istart + n, :)
849 : enddo
850 : endif
851 316860 : istart = istart + nindxm(l, itype)
852 : END DO
853 : END DO
854 : endif
855 : END DO
856 : END DO
857 : ! call timestop("desymmetrize")
858 7326 : END SUBROUTINE desymmetrize
859 :
860 : ! bra_trafo1 rotates cprod at kpts%bkp(ikpt)(<=> not irreducible k-point) to cprod at ikpt (bkp(kpts%bkp(ikpt))), which is the
861 : ! symmetrie equivalent one
862 : ! isym maps kpts%bkp(ikpt) on ikpt
863 :
864 16 : subroutine bra_trafo(fi, mpdata, hybdat, nbands, ikpt, psize, phase, vecin, vecout)
865 : use m_types_fleurinput
866 : USE m_types_mpdata
867 : USE m_types_hybdat
868 : USE m_types_mat
869 : use m_constants
870 : use m_judft
871 : implicit none
872 : type(t_fleurinput), intent(in) :: fi
873 : type(t_mpdata), intent(in) :: mpdata
874 : TYPE(t_hybdat), INTENT(IN) :: hybdat
875 : INTEGER, INTENT(IN) :: ikpt, nbands, psize
876 : type(t_mat), INTENT(IN) :: vecin
877 : type(t_mat), INTENT(INOUT) :: vecout
878 : COMPLEX, INTENT(INOUT) :: phase(:, :)
879 :
880 16 : if (vecin%l_real) then
881 12 : call bra_trafo_real(fi, mpdata, hybdat, nbands, ikpt, psize, phase, vecin%data_r, vecout%data_r)
882 : else
883 4 : phase = cmplx_1
884 4 : call bra_trafo_cmplx(fi, mpdata, hybdat, nbands, ikpt, psize, vecin%data_c, vecout%data_c)
885 : end if
886 :
887 16 : end subroutine bra_trafo
888 :
889 12 : subroutine bra_trafo_real(fi, mpdata, hybdat, nbands, ikpt, psize, phase, matin_r, matout_r)
890 : use m_types_fleurinput
891 : USE m_types_mpdata
892 : USE m_types_hybdat
893 : use m_constants
894 : use m_judft
895 : implicit none
896 : type(t_fleurinput), intent(in) :: fi
897 : type(t_mpdata), intent(in) :: mpdata
898 : TYPE(t_hybdat), INTENT(IN) :: hybdat
899 : INTEGER, INTENT(IN) :: ikpt, nbands, psize
900 : REAL, INTENT(IN) :: matin_r(:, :)
901 : REAL, INTENT(INOUT) :: matout_r(:, :)
902 : COMPLEX, INTENT(INOUT) :: phase(:, :)
903 :
904 12 : COMPLEX, ALLOCATABLE :: vecin(:, :), vecout(:, :)
905 : integer :: ok, i, j, cnt
906 12 : integer :: igptm2_list(mpdata%n_g(ikpt))
907 :
908 6136 : phase = cmplx_0
909 12 : call timestart("bra trafo real")
910 :
911 28 : IF (maxval(fi%hybinp%lcutm1) > fi%atoms%lmaxd) call judft_error('bra_trafo: maxlcutm > atoms%lmaxd') ! very improbable case
912 12 : call find_corresponding_g(fi%sym, fi%kpts, mpdata, ikpt, igptm2_list)
913 :
914 : ! transform back to unsymmetrized product basis in case of inversion symmetry
915 : !$OMP parallel default(none) private(i,j, cnt, vecin, vecout, ok) &
916 12 : !$OMP shared(nbands, psize, fi, hybdat, mpdata, phase, matin_r, matout_r, ikpt, igptm2_list)
917 : allocate (vecin(size(matin_r, dim=1), 1), vecout(size(matin_r, dim=1), 1), stat=ok, source=cmplx_0)
918 : IF (ok /= 0) call judft_error('bra_trafo: error allocating vecin or vecout')
919 :
920 : !$OMP do collapse(2)
921 : DO i = 1, nbands
922 : DO j = 1, psize
923 : cnt = (i-1) * psize + j
924 : vecin(:,1) = matin_r(:,cnt)
925 : CALL desymmetrize(vecin(:hybdat%nbasp, 1), hybdat%nbasp, 1, &
926 : fi%atoms, fi%hybinp%lcutm1, maxval(fi%hybinp%lcutm1), mpdata%num_radbasfn, fi%sym)
927 :
928 : call bra_trafo_core(1, ikpt, 1, fi%sym, mpdata, &
929 : fi%hybinp, hybdat, fi%kpts, fi%atoms, igptm2_list, vecin(:,1:1), vecout(:,1:1))
930 :
931 : CALL symmetrize(vecout(:, 1:1), hybdat%nbasm(ikpt), 1, 1, &
932 : fi%atoms, fi%hybinp%lcutm1, maxval(fi%hybinp%lcutm1), mpdata%num_radbasfn, fi%sym)
933 :
934 : phase(j, i) = commonphase(vecout(:, 1), hybdat%nbasm(ikpt))
935 : matout_r(:, cnt) = real(vecout(:, 1)/phase(j, i))
936 : IF (any(abs(aimag(vecout(:, 1)/phase(j, i))) > 1e-8)) THEN
937 : WRITE (*, *) vecout(:, 1)/phase(j, i)
938 : call judft_error('bra_trafo: Residual imaginary part.')
939 : END IF
940 :
941 : END DO
942 : END DO
943 : !$OMP end do
944 :
945 : deallocate (vecout, vecin)
946 : !$OMP end parallel
947 :
948 12 : call timestop("bra trafo real")
949 12 : end subroutine bra_trafo_real
950 :
951 4 : subroutine bra_trafo_cmplx(fi, mpdata, hybdat, nbands, ikpt, psize, vecin_c, vecout_c)
952 : use m_constants
953 : use m_judft
954 : use m_types_fleurinput
955 : USE m_types_mpdata
956 : USE m_types_hybdat
957 : implicit none
958 : type(t_fleurinput), intent(in) :: fi
959 : type(t_mpdata), intent(in) :: mpdata
960 : TYPE(t_hybdat), INTENT(IN) :: hybdat
961 : INTEGER, INTENT(IN) :: ikpt, nbands, psize
962 : COMPLEX, INTENT(IN) :: vecin_c(:, :)
963 : COMPLEX, INTENT(INOUT) :: vecout_c(:, :)
964 :
965 4 : integer :: igptm2_list(mpdata%n_g(ikpt))
966 :
967 4 : call timestart("bra trafo cmplx")
968 :
969 12 : IF (maxval(fi%hybinp%lcutm1) > fi%atoms%lmaxd) call judft_error('bra_trafo: maxlcutm > fi%atoms%lmaxd') ! very improbable case
970 4 : call find_corresponding_g(fi%sym, fi%kpts, mpdata, ikpt, igptm2_list)
971 :
972 4 : call bra_trafo_core(nbands, ikpt, psize, fi%sym, mpdata, fi%hybinp, hybdat, fi%kpts, fi%atoms, igptm2_list, vecin_c, vecout_c)
973 :
974 4 : call timestop("bra trafo cmplx")
975 4 : end subroutine bra_trafo_cmplx
976 :
977 4910 : subroutine bra_trafo_core(nbands, ikpt, psize, sym, &
978 4910 : mpdata, hybinp, hybdat, kpts, atoms, igptm2_list, vecin1, vecout1)
979 : use m_constants
980 : USE m_types_mpdata
981 : USE m_types_hybinp
982 : use m_types_hybdat
983 : USE m_types_sym
984 : USE m_types_kpts
985 : USE m_types_atoms
986 : implicit none
987 : type(t_mpdata), intent(in) :: mpdata
988 : TYPE(t_hybinp), INTENT(IN) :: hybinp
989 : TYPE(t_hybdat), INTENT(IN) :: hybdat
990 : TYPE(t_sym), INTENT(IN) :: sym
991 : TYPE(t_kpts), INTENT(IN) :: kpts
992 : TYPE(t_atoms), INTENT(IN) :: atoms
993 : integer, intent(in) :: igptm2_list(:)
994 :
995 : INTEGER, INTENT(IN) :: ikpt, nbands, psize
996 :
997 : COMPLEX, intent(in) :: vecin1(:, :)
998 : complex, intent(inout) :: vecout1(:, :)
999 :
1000 : INTEGER :: nrkpt, itype, ic, l, n, i, nn, i1, i2, j1, j2
1001 : INTEGER :: igptm, igptm2, igptp, iiop, inviop
1002 : COMPLEX :: cexp, cdum
1003 :
1004 : INTEGER :: rrot(3, 3), invrot(3, 3)
1005 48044 : INTEGER :: pnt(maxval(mpdata%num_radbasfn), 0:maxval(hybinp%lcutm1), atoms%nat)
1006 : INTEGER :: g(3), g1(3)
1007 : REAL :: rkpt(3), rkpthlp(3), trans(3)
1008 : COMPLEX :: dwgn(-maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), &
1009 35720 : -maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), 0:maxval(hybinp%lcutm1))
1010 :
1011 : ! call timestart("bra_trafo_core")
1012 : ! call timestart("setup")
1013 4910 : IF (kpts%bksym(ikpt) <= sym%nop) THEN
1014 4910 : inviop = sym%invtab(kpts%bksym(ikpt))
1015 63830 : rrot = transpose(sym%mrot(:, :, sym%invtab(kpts%bksym(ikpt))))
1016 4910 : invrot = sym%mrot(:, :, sym%invtab(kpts%bksym(ikpt)))
1017 19640 : trans = sym%tau(:, kpts%bksym(ikpt))
1018 :
1019 : dwgn(-maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), -maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), 0:maxval(hybinp%lcutm1)) &
1020 2284500 : = hybinp%d_wgn2(-maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), -maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), 0:maxval(hybinp%lcutm1), inviop)
1021 :
1022 : ELSE
1023 0 : iiop = kpts%bksym(ikpt) - sym%nop
1024 0 : inviop = sym%invtab(iiop) + sym%nop
1025 0 : rrot = -transpose(sym%mrot(:, :, sym%invtab(iiop)))
1026 : invrot = sym%mrot(:, :, sym%invtab(iiop))
1027 0 : trans = sym%tau(:, iiop)
1028 :
1029 : dwgn(-maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), -maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), 0:maxval(hybinp%lcutm1)) &
1030 0 : = conjg(hybinp%d_wgn2(-maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), -maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), 0:maxval(hybinp%lcutm1), inviop))
1031 : END IF
1032 :
1033 122750 : rkpt = matmul(rrot, kpts%bkf(:, kpts%bkp(ikpt)))
1034 4910 : rkpthlp = rkpt
1035 19640 : rkpt = kpts%to_first_bz(rkpt)
1036 19640 : g = nint(rkpthlp - rkpt)
1037 : ! call timestop("setup")
1038 :
1039 : !test
1040 : ! call timestart("test")
1041 4910 : nrkpt = 0
1042 27346 : DO i = 1, kpts%nkptf
1043 109384 : IF (maxval(abs(rkpt - kpts%bkf(:, i))) <= 1E-06) THEN
1044 : nrkpt = i
1045 : EXIT
1046 : END IF
1047 : END DO
1048 4910 : IF (nrkpt /= ikpt) THEN
1049 0 : PRINT *, kpts%bkp(ikpt), ikpt
1050 0 : PRINT *, kpts%bkf(:, ikpt)
1051 0 : PRINT *, kpts%bkf(:, kpts%bkp(ikpt))
1052 0 : PRINT *, rkpt
1053 :
1054 0 : call judft_error('bra_trafo: rotation failed')
1055 : END IF
1056 : ! call timestop("test")
1057 :
1058 : ! Define pointer to first mixed-basis functions (with m = -l)
1059 : ! call timestart("def pointer to first mpb")
1060 : i = 0
1061 11072 : do ic = 1, atoms%nat
1062 6162 : itype = atoms%itype(ic)
1063 41882 : DO l = 0, hybinp%lcutm1(itype)
1064 269708 : DO n = 1, mpdata%num_radbasfn(l, itype)
1065 238898 : i = i + 1
1066 269708 : pnt(n, l, ic) = i
1067 : END DO
1068 36972 : i = i + mpdata%num_radbasfn(l, itype)*2*l
1069 : END DO
1070 : END DO
1071 : ! call timestop("def pointer to first mpb")
1072 :
1073 : ! Multiplication
1074 : ! MT
1075 : ! call timestart("MT part")
1076 19640 : cexp = exp(ImagUnit*tpi_const*dot_product(kpts%bkf(:, ikpt) + g, trans(:)))
1077 : !$OMP parallel do default(none) private(ic, itype, cdum, l, nn, n, i1, i2, j1, j2, i)&
1078 4910 : !$OMP shared(atoms, cexp, hybinp, kpts, mpdata, pnt, dwgn, vecin1, vecout1, ikpt, g, nbands, psize)
1079 : do ic = 1, atoms%nat
1080 : itype = atoms%itype(ic)
1081 :
1082 : cdum = cexp*exp(-ImagUnit*tpi_const*dot_product(g, atoms%taual(:, hybinp%map(ic, kpts%bksym(ikpt)))))
1083 :
1084 : DO l = 0, hybinp%lcutm1(itype)
1085 : nn = mpdata%num_radbasfn(l, itype)
1086 : DO n = 1, nn
1087 :
1088 : i1 = pnt(n, l, ic)
1089 : i2 = i1 + nn*2*l
1090 : j1 = pnt(n, l, hybinp%map(ic, kpts%bksym(ikpt)))
1091 : j2 = j1 + nn*2*l
1092 :
1093 : DO i = 1, nbands*psize
1094 : vecout1(i1:i2:nn, i) = cdum*matmul(vecin1(j1:j2:nn, i), dwgn(-l:l, -l:l, l))
1095 : END DO
1096 : END DO
1097 : END DO
1098 : END DO
1099 : !$OMP end parallel do
1100 : ! call timestop("MT part")
1101 :
1102 : ! PW
1103 : ! call timestart("PW part")
1104 : !$OMP parallel do default(none) private(igptm, igptp, g1, igptm2, i, cdum) &
1105 4910 : !$OMP shared(vecout1, vecin1, mpdata, ikpt, igptm2_list, kpts, rrot, g, hybdat, trans, nbands, psize)
1106 : DO igptm = 1, mpdata%n_g(kpts%bkp(ikpt))
1107 : igptp = mpdata%gptm_ptr(igptm, kpts%bkp(ikpt))
1108 : g1 = matmul(rrot, mpdata%g(:, igptp)) + g
1109 : igptm2 = igptm2_list(igptm)
1110 :
1111 : cdum = exp(ImagUnit*tpi_const*dot_product(kpts%bkf(:, ikpt) + g1, trans(:)))
1112 : vecout1(hybdat%nbasp + igptm, :) = cdum*vecin1(hybdat%nbasp + igptm2, :)
1113 : END DO
1114 : !$OMP end parallel do
1115 : ! call timestop("PW part")
1116 : ! call timestop("bra_trafo_core")
1117 4910 : end subroutine bra_trafo_core
1118 :
1119 16 : subroutine find_corresponding_g(sym, kpts, mpdata, ikpt, igptm2_list)
1120 : use m_types_sym
1121 : USE m_types_kpts
1122 : USE m_types_mpdata
1123 : implicit none
1124 : type(t_sym), intent(in) :: sym
1125 : type(t_kpts), intent(in) :: kpts
1126 : type(t_mpdata), intent(in) :: mpdata
1127 : integer, intent(in) :: ikpt
1128 : integer, intent(inout) :: igptm2_list(:)
1129 :
1130 : integer :: igptm, igptp, g1(3), igptm2, i, iiop
1131 : integer :: g(3), rrot(3, 3)
1132 : REAL :: rkpt(3), rkpthlp(3)
1133 :
1134 16 : call timestart("find corresponding g")
1135 16 : call timestart("setup")
1136 16 : IF (kpts%bksym(ikpt) <= sym%nop) THEN
1137 208 : rrot = transpose(sym%mrot(:, :, sym%invtab(kpts%bksym(ikpt))))
1138 : ELSE
1139 0 : iiop = kpts%bksym(ikpt) - sym%nop
1140 0 : rrot = -transpose(sym%mrot(:, :, sym%invtab(iiop)))
1141 : END IF
1142 :
1143 400 : rkpt = matmul(rrot, kpts%bkf(:, kpts%bkp(ikpt)))
1144 16 : rkpthlp = rkpt
1145 64 : rkpt = kpts%to_first_bz(rkpt)
1146 64 : g = nint(rkpthlp - rkpt)
1147 16 : call timestop("setup")
1148 :
1149 : !$OMP parallel do default(none) schedule(dynamic, 10) private(igptm, igptp, g1, igptm2) &
1150 16 : !$OMP shared(kpts, mpdata, ikpt, rrot, g, igptm2_list)
1151 : do igptm = 1, mpdata%n_g(kpts%bkp(ikpt))
1152 : igptp = mpdata%gptm_ptr(igptm, kpts%bkp(ikpt))
1153 : g1 = matmul(rrot, mpdata%g(:, igptp)) + g
1154 :
1155 : igptm2 = 0
1156 : DO i = 1, mpdata%n_g(ikpt)
1157 : IF (maxval(abs(g1 - mpdata%g(:, mpdata%gptm_ptr(i, ikpt)))) <= 1E-06) THEN
1158 : igptm2 = i
1159 : EXIT
1160 : END IF
1161 : END DO
1162 : IF (igptm2 == 0) THEN
1163 : WRITE (*, *) kpts%bkp(ikpt), ikpt, g1
1164 : WRITE (*, *) mpdata%n_g(kpts%bkp(ikpt)), mpdata%n_g(ikpt)
1165 : WRITE (*, *)
1166 : WRITE (*, *) igptp, mpdata%g(:, igptp)
1167 : WRITE (*, *) g
1168 : WRITE (*, *) rrot
1169 : WRITE (*, *) "Failed tests:", g1
1170 : DO i = 1, mpdata%n_g(ikpt)
1171 : WRITE (*, *) mpdata%g(:, mpdata%gptm_ptr(i, ikpt))
1172 : END DO
1173 : call judft_error('bra_trafo: G-point not found in G-point set.')
1174 : END IF
1175 :
1176 : igptm2_list(igptm) = igptm2
1177 : enddo
1178 : !$OMP end parallel do
1179 16 : call timestop("find corresponding g")
1180 16 : end subroutine find_corresponding_g
1181 :
1182 : ! Determines common phase factor (with unit norm)
1183 5066 : function commonphase(carr, n) result(cfac)
1184 : USE m_juDFT
1185 : IMPLICIT NONE
1186 : INTEGER, INTENT(IN) :: n
1187 : COMPLEX, INTENT(IN) :: carr(n)
1188 : COMPLEX :: cfac
1189 : REAL :: rdum, rmax
1190 : INTEGER :: i
1191 :
1192 5066 : cfac = 0
1193 5066 : rmax = 0
1194 49083 : DO i = 1, n
1195 49083 : rdum = abs(carr(i))
1196 49083 : IF (rdum > 1e-6) THEN
1197 5066 : cfac = carr(i)/rdum
1198 5066 : EXIT
1199 44017 : ELSE IF (rdum > rmax) THEN
1200 17390 : cfac = carr(i)/rdum
1201 17390 : rmax = rdum
1202 : END IF
1203 : END DO
1204 33658 : IF (abs(cfac) < 1e-10 .and. all(abs(carr) > 1e-10)) THEN
1205 0 : WRITE (999, *) carr
1206 0 : call judft_error('commonphase: Could not determine common phase factor. (Wrote carr to fort.999)')
1207 : END IF
1208 5066 : END function commonphase
1209 :
1210 0 : function commonphase_mtx(mtx, dim1, dim2) result(cfac)
1211 : implicit none
1212 :
1213 : COMPLEX, INTENT(IN) :: mtx(:, :)
1214 : integer, intent(in) :: dim1, dim2
1215 : COMPLEX :: cfac
1216 : REAL :: rdum, rmax
1217 : INTEGER :: i, j
1218 :
1219 0 : do j = 1, dim2
1220 0 : do i = 1, dim1
1221 0 : rdum = abs(mtx(i, j))
1222 0 : IF (rdum > 1e-6) THEN
1223 0 : cfac = mtx(i, j)/rdum
1224 0 : EXIT
1225 0 : ELSE IF (rdum > rmax) THEN
1226 0 : cfac = mtx(i, j)/rdum
1227 0 : rmax = rdum
1228 : END IF
1229 : end do
1230 : end do
1231 :
1232 0 : IF (abs(cfac) < 1e-10 .and. all(abs(mtx(:dim1, :dim2)) > 1e-10)) THEN
1233 0 : WRITE (999, *) mtx(:dim1, :dim2)
1234 0 : call judft_error('commonphase: Could not determine common phase factor. (Wrote carr to fort.999)')
1235 : END IF
1236 0 : END function commonphase_mtx
1237 :
1238 14036916 : SUBROUTINE bramat_trafo(vecin, igptm_in, ikpt0, iop, writevec, pointer, sym, &
1239 11556 : rrot, invrrot, mpdata, hybinp, kpts, maxlcutm, atoms, lcutm, nindxm, maxindxm, &
1240 104004 : dwgn, nbasp, nbasm, vecout, igptm_out)
1241 :
1242 : USE m_constants
1243 : USE m_util
1244 : USE m_types_mpdata
1245 : USE m_types_hybinp
1246 : USE m_types_sym
1247 : USE m_types_kpts
1248 : USE m_types_atoms
1249 : IMPLICIT NONE
1250 : type(t_mpdata), intent(in) :: mpdata
1251 : TYPE(t_hybinp), INTENT(IN) :: hybinp
1252 : TYPE(t_sym), INTENT(IN) :: sym
1253 : TYPE(t_kpts), INTENT(IN) :: kpts
1254 : TYPE(t_atoms), INTENT(IN) :: atoms
1255 :
1256 : ! - scalars
1257 : INTEGER, INTENT(IN) :: ikpt0, igptm_in, iop, maxindxm
1258 : INTEGER, INTENT(IN) :: maxlcutm
1259 : INTEGER, INTENT(IN) :: nbasp
1260 : LOGICAL, INTENT(IN) :: writevec
1261 : INTEGER, INTENT(INOUT) :: igptm_out
1262 : ! - arrays -
1263 : INTEGER, INTENT(IN) :: rrot(:, :), invrrot(:, :)
1264 : INTEGER, INTENT(IN) :: lcutm(atoms%ntype), &
1265 : nindxm(0:maxlcutm, atoms%ntype)
1266 : INTEGER, INTENT(IN) :: nbasm(:)
1267 : INTEGER, INTENT(IN) :: pointer( &
1268 : minval(mpdata%g(1, :)) - 1:maxval(mpdata%g(1, :)) + 1, &
1269 : minval(mpdata%g(2, :)) - 1:maxval(mpdata%g(2, :)) + 1, &
1270 : minval(mpdata%g(3, :)) - 1:maxval(mpdata%g(3, :)) + 1)
1271 :
1272 : COMPLEX, INTENT(IN) :: vecin(:)
1273 : COMPLEX, INTENT(IN) :: dwgn(-maxlcutm:maxlcutm, &
1274 : -maxlcutm:maxlcutm, &
1275 : 0:maxlcutm)
1276 : COMPLEX, INTENT(INOUT) :: vecout(maxval(nbasm), 1)
1277 :
1278 : ! - private scalars -
1279 : INTEGER :: itype, ieq, ic, l, n, i, nn, i1, i2, j1, j2
1280 : INTEGER :: igptm, igptm2, igptp, isym
1281 : INTEGER :: ikpt1
1282 : LOGICAL :: trs, touch
1283 : COMPLEX :: cexp, cdum
1284 : ! - private arrays -
1285 11556 : INTEGER :: pnt(maxindxm, 0:maxlcutm, atoms%nat), g(3), &
1286 104004 : g1(3), iarr(maxval(mpdata%n_g))
1287 : REAL :: rkpt(3), rkpthlp(3), trans(3)
1288 11556 : COMPLEX :: vecin1(nbasm(ikpt0))
1289 104004 : COMPLEX :: carr(maxval(mpdata%n_g))
1290 :
1291 11556 : call timestart("bramat_trafo")
1292 5751064 : igptm_out = -1; vecout = CMPLX_NOT_INITALIZED; touch = .false.
1293 :
1294 11556 : IF (iop <= sym%nop) THEN
1295 10632 : isym = iop
1296 10632 : trs = .false.
1297 42528 : trans = sym%tau(:, isym)
1298 : ELSE
1299 924 : isym = iop - sym%nop
1300 924 : trs = .true.
1301 3696 : trans = sym%tau(:, isym)
1302 : END IF
1303 :
1304 288900 : rkpthlp = matmul(rrot, kpts%bkf(:, ikpt0))
1305 11556 : rkpt = kpts%to_first_bz(rkpthlp)
1306 46224 : g = nint(rkpthlp - rkpt)
1307 : !
1308 : ! determine number of rotated k-point bk(:,ikpt) -> ikpt1
1309 : !
1310 11556 : call timestart("det. kpoint")
1311 20108 : DO i = 1, kpts%nkpt
1312 80432 : IF (maxval(abs(rkpt - kpts%bkf(:, i))) <= 1E-06) THEN
1313 : ikpt1 = i
1314 : EXIT
1315 : END IF
1316 : END DO
1317 11556 : call timestop("det. kpoint")
1318 :
1319 11556 : call timestart("calc igptm_out")
1320 896060 : DO igptm = 1, mpdata%n_g(ikpt1)
1321 896060 : igptp = mpdata%gptm_ptr(igptm, ikpt1)
1322 14336960 : g1 = matmul(invrrot, mpdata%g(:, igptp) - g)
1323 896060 : igptm2 = pointer(g1(1), g1(2), g1(3))
1324 896060 : IF (igptm2 == igptm_in) THEN
1325 11556 : igptm_out = igptm
1326 11556 : touch = .true.
1327 11556 : IF (writevec) THEN
1328 13936 : cdum = exp(ImagUnit*tpi_const*dot_product(kpts%bkf(:, ikpt1) + mpdata%g(:, igptp), trans))
1329 : EXIT
1330 : ELSE
1331 8072 : call timestop("calc igptm_out")
1332 8072 : call timestop("bramat_trafo")
1333 : RETURN
1334 : END IF
1335 : END IF
1336 : END DO
1337 3484 : call timestop("calc igptm_out")
1338 :
1339 3484 : if (.not. touch) call judft_error("g-point could not be found.")
1340 : ! Transform back to unsymmetrized product basis in case of inversion symmetry.
1341 :
1342 3484 : call timestart("desymm")
1343 1713832 : vecout(:nbasm(ikpt0), 1) = vecin(:nbasm(ikpt0))
1344 3484 : if (sym%invs) CALL desymmetrize(vecout, nbasp, 1, &
1345 2420 : atoms, lcutm, maxlcutm, nindxm, sym)
1346 3484 : call timestop("desymm")
1347 :
1348 : ! Right-multiplication
1349 : ! PW
1350 3484 : call timestart("right-multiply")
1351 174232 : IF (trs) THEN; vecin1(:nbasm(ikpt0)) = cdum*conjg(vecout(:nbasm(ikpt0), 1))
1352 1542760 : ELSE; vecin1(:nbasm(ikpt0)) = cdum*vecout(:nbasm(ikpt0), 1)
1353 : END IF
1354 3484 : call timestop("right-multiply")
1355 :
1356 : ! Define pointer to first mixed-basis functions (with m = -l)
1357 3484 : call timestart("def pointer")
1358 3484 : i = 0
1359 3484 : ic = 0
1360 10020 : DO itype = 1, atoms%ntype
1361 16556 : DO ieq = 1, atoms%neq(itype)
1362 6536 : ic = ic + 1
1363 45752 : DO l = 0, lcutm(itype)
1364 274532 : DO n = 1, nindxm(l, itype)
1365 241852 : i = i + 1
1366 274532 : pnt(n, l, ic) = i
1367 : END DO
1368 39216 : i = i + nindxm(l, itype)*2*l
1369 : END DO
1370 : END DO
1371 : END DO
1372 3484 : call timestop("def pointer")
1373 :
1374 : ! Left-multiplication
1375 : ! MT
1376 3484 : call timestart("left multi MT")
1377 13936 : cexp = exp(-ImagUnit*tpi_const*dot_product(kpts%bkf(:, ikpt1) + g, trans))
1378 3484 : ic = 0
1379 10020 : DO itype = 1, atoms%ntype
1380 16556 : DO ieq = 1, atoms%neq(itype)
1381 6536 : ic = ic + 1
1382 26144 : cdum = cexp*exp(ImagUnit*tpi_const*dot_product(g, atoms%taual(:, ic)))
1383 6536 : cdum = conjg(cdum)
1384 45752 : DO l = 0, lcutm(itype)
1385 32680 : nn = nindxm(l, itype)
1386 281068 : DO n = 1, nn
1387 :
1388 241852 : i1 = pnt(n, l, ic)
1389 241852 : i2 = i1 + nn*2*l
1390 241852 : j1 = pnt(n, l, hybinp%map(ic, sym%invtab(isym)))
1391 241852 : j2 = j1 + nn*2*l
1392 :
1393 10568468 : vecout(i1:i2:nn, 1) = cdum*matmul(dwgn(-l:l, -l:l, l), vecin1(j1:j2:nn))
1394 :
1395 : END DO
1396 : END DO
1397 : END DO
1398 : END DO
1399 3484 : call timestop("left multi MT")
1400 :
1401 : ! PW
1402 3484 : call timestart("left multi pw")
1403 602924 : DO igptm = 1, mpdata%n_g(ikpt1)
1404 599440 : igptp = mpdata%gptm_ptr(igptm, ikpt1)
1405 9591040 : g1 = matmul(invrrot, mpdata%g(:, igptp) - g)
1406 599440 : iarr(igptm) = pointer(g1(1), g1(2), g1(3))
1407 2401244 : carr(igptm) = exp(-ImagUnit*tpi_const*dot_product(kpts%bkf(:, ikpt1) + mpdata%g(:, igptp), trans))
1408 : END DO
1409 602924 : DO i1 = 1, mpdata%n_g(ikpt1)
1410 602924 : vecout(nbasp + i1, 1) = carr(i1)*vecin1(nbasp + iarr(i1))
1411 : END DO
1412 3484 : call timestop("left multi pw")
1413 :
1414 : ! If inversion symmetry is applicable, symmetrize to make the values real.
1415 3484 : call timestart("symmetrize")
1416 3484 : if (sym%invs) CALL symmetrize(vecout, nbasp, 1, 1, &
1417 2420 : atoms, lcutm, maxlcutm, nindxm, sym)
1418 3484 : call timestop("symmetrize")
1419 3484 : call timestop("bramat_trafo")
1420 : END SUBROUTINE bramat_trafo
1421 :
1422 1511224 : END MODULE m_trafo
|
