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_symmetrizeh
8 :
9 : ! symmetrize the Hamiltonian according to the symmetry operations of the little group of k
10 :
11 : CONTAINS
12 :
13 24 : SUBROUTINE symmetrizeh(atoms, bk, jsp, lapw, sym, cell, nsymop, psym, hmat)
14 : use m_map_to_unit
15 : USE m_juDFT
16 : USE m_types
17 : USE m_constants
18 :
19 : IMPLICIT NONE
20 :
21 : TYPE(t_sym), INTENT(IN) :: sym
22 : TYPE(t_cell), INTENT(IN) :: cell
23 : TYPE(t_atoms), INTENT(IN) :: atoms
24 : TYPE(t_lapw), INTENT(IN) :: lapw
25 : TYPE(t_mat), INTENT(INOUT) :: hmat
26 :
27 : ! scalars
28 : INTEGER, INTENT(IN) :: nsymop, jsp
29 :
30 : ! arrays
31 : INTEGER, INTENT(IN) :: psym(:)
32 : REAL, INTENT(IN) :: bk(:)
33 :
34 : ! local scalars
35 : INTEGER :: ilotot, itype, itype1, ilo, ilo1
36 : INTEGER :: iatom, iatom1, iiatom
37 : INTEGER :: i, ieq, ieq1, m
38 : INTEGER :: igpt_lo, igpt_lo1, igpt_lo2, igpt1_lo1
39 : INTEGER :: igpt1_lo2, isym, iop, ic, ic1, ic2
40 : INTEGER :: igpt, igpt1, igpt2, igpt3
41 : INTEGER :: invsfct, invsfct1, idum
42 : INTEGER :: l, l1, lm, j, ok, ratom, ratom1, nrgpt
43 : COMPLEX, PARAMETER :: img = (0.0, 1.0)
44 : COMPLEX :: cdum, cdum2
45 :
46 : ! local arrays
47 24 : INTEGER :: l_lo(atoms%nlotot)
48 24 : INTEGER :: itype_lo(atoms%nlotot)
49 24 : INTEGER :: gpt_lo(3, atoms%nlotot), gpthlp(3), g(3)
50 24 : INTEGER :: lo_indx(atoms%nlod, atoms%nat)
51 24 : INTEGER :: rot(3, 3, nsymop), rrot(3, 3, nsymop)
52 :
53 24 : INTEGER, ALLOCATABLE :: pointer_apw(:, :)
54 24 : INTEGER, ALLOCATABLE :: ipiv(:)
55 24 : INTEGER, ALLOCATABLE :: map(:, :)
56 :
57 : REAL :: rtaual(3), kghlp(3)
58 : REAL :: rotkpt(3)
59 24 : REAL :: trans(3, nsymop)
60 24 : COMPLEX, ALLOCATABLE :: c_lo(:, :, :, :), c_rot(:, :, :, :, :), y(:)
61 24 : COMPLEX, ALLOCATABLE :: cfac(:, :), chelp(:, :)
62 :
63 24 : LOGICAL :: ldum(lapw%nv(jsp) + atoms%nlotot, lapw%nv(jsp) + atoms%nlotot)
64 :
65 24 : TYPE(t_mat) :: hmatTemp
66 :
67 24 : CALL hmatTemp%init(hmat%l_real, hmat%matsize1, hmat%matsize2)
68 24 : IF (hmat%l_real) THEN
69 461460 : hmatTemp%data_r = hmat%data_r
70 : ELSE
71 196256 : hmatTemp%data_c = CONJG(hmat%data_c)
72 : END IF
73 :
74 : ! calculate rotations in reciprocal space
75 744 : DO isym = 1, nsymop
76 720 : iop = psym(isym)
77 744 : IF (iop <= sym%nop) THEN
78 8372 : rrot(:, :, isym) = TRANSPOSE(sym%mrot(:, :, sym%invtab(iop)))
79 : ELSE
80 988 : rrot(:, :, isym) = -TRANSPOSE(sym%mrot(:, :, sym%invtab(iop - sym%nop)))
81 : END IF
82 : END DO
83 :
84 : ! calculate rotations in real space (internal coordinates)
85 744 : DO isym = 1, nsymop
86 720 : iop = psym(isym)
87 744 : IF (iop <= sym%nop) THEN
88 8372 : rot(:, :, isym) = sym%mrot(:, :, iop)
89 2576 : trans(:, isym) = sym%tau(:, iop)
90 : ELSE
91 988 : rot(:, :, isym) = sym%mrot(:, :, iop - sym%nop)
92 304 : trans(:, isym) = sym%tau(:, iop - sym%nop)
93 : END IF
94 : END DO
95 :
96 : ! caclulate mapping of atoms
97 96 : allocate(map(nsymop, atoms%nat))
98 1084 : map = 0
99 24 : iatom = 0
100 60 : DO itype = 1, atoms%ntype
101 36 : iiatom = atoms%firstAtom(itype) - 1
102 96 : DO ieq = 1, atoms%neq(itype)
103 36 : iatom = iatom + 1
104 1096 : DO isym = 1, nsymop
105 28672 : rtaual = MATMUL(rot(:, :, isym), atoms%taual(:, iatom)) + trans(:, isym)
106 1024 : iatom1 = 0
107 2048 : DO ieq1 = 1, atoms%neq(itype)
108 9216 : IF (norm2(map_to_unit(rtaual - atoms%taual(:, iiatom + ieq1))) < 1e-6) THEN
109 3072 : iatom1 = iiatom + ieq1
110 : END IF
111 : END DO
112 1024 : IF (iatom1 == 0) call judft_error('symmetrizeh_new: error finding rotated atomic position')
113 1060 : map(isym, iatom) = iatom1
114 : END DO
115 : END DO
116 : END DO
117 :
118 : ! initialze pointer_apw and the apw part of cfac
119 168 : allocate(pointer_apw(lapw%nv(jsp), nsymop), cfac(lapw%nv(jsp) + atoms%nlotot, nsymop), stat=ok)
120 24 : IF (ok /= 0) call judft_error('symmetrizeh_new: failure allocation pointer_apw,cfac')
121 :
122 99864 : pointer_apw = 0
123 103656 : cfac = 0
124 :
125 744 : DO isym = 1, nsymop
126 :
127 : ! determine vector g, which map Rk back into BZ
128 : ! Rk - G = k => G = Rk-k
129 18000 : rotkpt = MATMUL(rrot(:, :, isym), bk(:))
130 2880 : g = NINT(rotkpt - bk)
131 :
132 99864 : DO igpt = 1, lapw%nv(jsp)
133 : !rotate G vector corresponding to isym
134 1585920 : gpthlp = MATMUL(rrot(:, :, isym), lapw%gvec(:, igpt, jsp)) + g
135 : ! determine number of gpthlp
136 99120 : nrgpt = 0
137 8516744 : DO i = 1, lapw%nv(jsp)
138 34066976 : IF (MAXVAL(ABS(gpthlp - lapw%gvec(:, i, jsp))) <= 1E-06) THEN
139 : nrgpt = i
140 : EXIT
141 : END IF
142 : END DO
143 99120 : IF (nrgpt == 0) THEN
144 0 : PRINT *, igpt
145 0 : PRINT *, lapw%gvec(:, igpt, jsp)
146 0 : PRINT *, gpthlp
147 0 : PRINT *, g
148 0 : PRINT *, bk
149 0 : DO i = 1, lapw%nv(jsp)
150 0 : WRITE (oUnit, *) i, lapw%gvec(:, i, jsp)
151 : ENDDO
152 0 : call judft_error('symmetrizeh_new: rotated G point not found')
153 : END IF
154 99120 : pointer_apw(igpt, isym) = nrgpt
155 397200 : cfac(igpt, isym) = EXP(-2*pi_const*img*(dot_PRODUCT(bk(:) + gpthlp(:), trans(:, isym))))
156 : END DO
157 : END DO
158 :
159 : ! average apw-part of symmetry-equivalent matrix elements
160 :
161 657692 : ldum = .TRUE.
162 3496 : DO i = 1, lapw%nv(jsp)
163 311800 : DO j = 1, i
164 : cdum = 0
165 : ic = 0
166 8825048 : DO isym = 1, nsymop
167 8516744 : iop = psym(isym)
168 8516744 : igpt = pointer_apw(i, isym)
169 8516744 : igpt1 = pointer_apw(j, isym)
170 :
171 8825048 : IF (iop <= sym%nop) THEN
172 7413272 : IF ((igpt /= 0) .AND. (igpt1 /= 0)) THEN
173 : ic = ic + 1
174 7413272 : IF (hmatTemp%l_real) THEN
175 6309800 : cdum = cdum + CONJG(cfac(i, isym))*hmatTemp%data_r(igpt1, igpt)*cfac(j, isym)
176 : ELSE
177 1103472 : cdum = cdum + CONJG(cfac(i, isym))*hmatTemp%data_c(igpt1, igpt)*cfac(j, isym)
178 : END IF
179 : END IF
180 : ELSE
181 1103472 : IF ((igpt /= 0) .AND. (igpt1 /= 0)) THEN
182 : ic = ic + 1
183 1103472 : IF (hmatTemp%l_real) THEN
184 0 : cdum = cdum + CONJG(CONJG(cfac(i, isym))*hmatTemp%data_r(igpt1, igpt)*cfac(j, isym))
185 : ELSE
186 1103472 : cdum = cdum + CONJG(CONJG(cfac(i, isym))*hmatTemp%data_c(igpt1, igpt)*cfac(j, isym))
187 : END IF
188 : END IF
189 : END IF
190 : END DO
191 :
192 8825048 : cdum = cdum!/ic
193 8828520 : DO isym = 1, nsymop
194 8516744 : iop = psym(isym)
195 8516744 : igpt = pointer_apw(i, isym)
196 8516744 : igpt1 = pointer_apw(j, isym)
197 8516744 : IF ((igpt == 0) .OR. (igpt1 == 0)) CYCLE
198 8516744 : IF (igpt1 > igpt) CYCLE
199 8328860 : IF (ldum(igpt, igpt1)) THEN
200 308304 : IF (hmat%l_real) THEN
201 220732 : IF (iop <= sym%nop) THEN
202 220732 : hmat%data_r(igpt1, igpt) = real(cdum/(CONJG(cfac(i, isym))*cfac(j, isym)))
203 220732 : ldum(igpt, igpt1) = .FALSE.
204 : ELSE
205 0 : hmat%data_r(igpt1, igpt) = real(CONJG(cdum/(CONJG(cfac(i, isym))*cfac(j, isym))))
206 0 : ldum(igpt, igpt1) = .FALSE.
207 : END IF
208 220732 : hmat%data_r(igpt, igpt1) = hmat%data_r(igpt1, igpt)
209 : ELSE
210 87572 : IF (iop <= sym%nop) THEN
211 46016 : hmat%data_c(igpt1, igpt) = cdum/(CONJG(cfac(i, isym))*cfac(j, isym))
212 46016 : ldum(igpt, igpt1) = .FALSE.
213 : ELSE
214 41556 : hmat%data_c(igpt1, igpt) = CONJG(cdum/(CONJG(cfac(i, isym))*cfac(j, isym)))
215 41556 : ldum(igpt, igpt1) = .FALSE.
216 : END IF
217 87572 : hmat%data_c(igpt, igpt1) = CONJG(hmat%data_c(igpt1, igpt))
218 : END IF
219 : END IF
220 : END DO
221 : END DO
222 : END DO
223 :
224 : ! average lo-part of matrix elements
225 48 : IF (ANY(atoms%nlo /= 0)) THEN
226 :
227 : ! preparations
228 : ilotot = 0
229 : iatom = 0
230 60 : DO itype = 1, atoms%ntype
231 96 : DO ieq = 1, atoms%neq(itype)
232 36 : iatom = iatom + 1
233 72 : IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN
234 36 : IF (sym%invsat(iatom) == 0) invsfct = 1
235 36 : IF (sym%invsat(iatom) == 1) invsfct = 2
236 84 : DO ilo = 1, atoms%nlo(itype)
237 48 : l = atoms%llo(ilo, itype)
238 216 : DO m = 1, invsfct*(2*l + 1)
239 132 : ilotot = ilotot + 1
240 132 : l_lo(ilotot) = l
241 132 : itype_lo(ilotot) = itype
242 576 : gpt_lo(:, ilotot) = lapw%gvec(:, lapw%kvec(m,ilo,iatom), jsp)
243 : END DO
244 : END DO
245 : END IF
246 : END DO
247 : END DO
248 :
249 : ! calculate expansion coefficients for local orbitals
250 24 : IF (hmat%l_real) THEN
251 180 : allocate(c_lo(4*MAXVAL(l_lo) + 2, 4*MAXVAL(l_lo) + 2, atoms%nlod, atoms%nat), stat=ok)
252 : ELSE
253 96 : allocate(c_lo(2*MAXVAL(l_lo) + 1, 2*MAXVAL(l_lo) + 1, atoms%nlod, atoms%nat), stat=ok)
254 : END IF
255 :
256 24 : IF (ok /= 0) call judft_error('symmetrizeh_new: failure allocation c_lo')
257 :
258 60 : iatom = 0
259 60 : ilotot = 0
260 120 : lo_indx = 0
261 60 : DO itype = 1, atoms%ntype
262 96 : DO ieq = 1, atoms%neq(itype)
263 36 : iatom = iatom + 1
264 72 : IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN
265 36 : IF (sym%invsat(iatom) == 0) invsfct = 1
266 36 : IF (sym%invsat(iatom) == 1) invsfct = 2
267 :
268 84 : DO ilo = 1, atoms%nlo(itype)
269 48 : l = atoms%llo(ilo, itype)
270 144 : allocate(y((l + 1)**2))
271 48 : lo_indx(ilo, iatom) = ilotot + 1
272 :
273 180 : DO igpt_lo = 1, invsfct*(2*l + 1)
274 132 : ilotot = ilotot + 1
275 528 : kghlp = bk(:) + gpt_lo(:, ilotot)
276 :
277 : !generate spherical harmonics
278 1716 : CALL harmonicsr(y, MATMUL(kghlp, cell%bmat), l)
279 :
280 132 : lm = l**2
281 528 : cdum = EXP(2*pi_const*img*dot_PRODUCT(kghlp, atoms%taual(:, iatom)))
282 :
283 180 : IF (invsfct == 1) THEN
284 612 : DO m = 1, 2*l + 1
285 480 : lm = lm + 1
286 612 : c_lo(m, igpt_lo, ilo, iatom) = cdum*CONJG(y(lm))
287 : END DO
288 : ELSE
289 0 : DO m = 1, 2*l + 1
290 0 : lm = lm + 1
291 0 : c_lo(m, igpt_lo, ilo, iatom) = cdum*CONJG(y(lm))
292 0 : c_lo(4*l + 3 - m, igpt_lo, ilo, iatom) = (-1)**(m - 1)*CONJG(cdum)*y(lm)
293 : END DO
294 : END IF
295 : END DO
296 84 : deallocate(y)
297 : END DO
298 : END IF
299 : END DO
300 : END DO
301 :
302 24 : IF (ilotot /= atoms%nlotot) call judft_error('symmetrizeh_new: failure counting local orbitals(ilotot)')
303 :
304 24 : IF (hmat%l_real) THEN
305 198 : allocate(c_rot(4*MAXVAL(l_lo) + 2, 4*MAXVAL(l_lo) + 2, atoms%nlod, atoms%nat, nsymop))
306 144 : allocate(chelp(4*MAXVAL(l_lo) + 2, 4*MAXVAL(l_lo) + 2))
307 : ELSE
308 102 : allocate(c_rot(2*MAXVAL(l_lo) + 1, 2*MAXVAL(l_lo) + 1, atoms%nlod, atoms%nat, nsymop))
309 84 : allocate(chelp(2*MAXVAL(l_lo) + 1, 2*MAXVAL(l_lo) + 1))
310 : END IF
311 :
312 744 : DO isym = 1, nsymop
313 : ! determine vector g, which map Rk back into BZ
314 : ! Rk - G = k => G = Rk-k
315 18000 : rotkpt = MATMUL(rrot(:, :, isym), bk(:))
316 2880 : g = NINT(rotkpt - bk)
317 :
318 : ilotot = 0
319 : iatom = 0
320 1768 : DO itype = 1, atoms%ntype
321 2768 : DO ieq = 1, atoms%neq(itype)
322 1024 : iatom = iatom + 1
323 1024 : ratom = map(isym, iatom)
324 :
325 2048 : IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN
326 1024 : IF (sym%invsat(iatom) == 0) invsfct = 1
327 1024 : IF (sym%invsat(iatom) == 1) THEN
328 0 : IF (sym%invsat(ratom) == 2) THEN
329 0 : ratom = sym%invsatnr(ratom)
330 : END IF
331 : invsfct = 2
332 : END IF
333 :
334 2464 : DO ilo = 1, atoms%nlo(itype)
335 1440 : l = atoms%llo(ilo, itype)
336 4320 : allocate(y((l + 1)**2))
337 :
338 5232 : DO igpt_lo = 1, invsfct*(2*l + 1)
339 3792 : ilotot = ilotot + 1
340 :
341 : !rotate G_lo corresponding to iop
342 60672 : gpthlp = MATMUL(rrot(:, :, isym), gpt_lo(:, ilotot)) + g
343 60672 : kghlp = bk(:) + MATMUL(rrot(:, :, isym), gpt_lo(:, ilotot)) + g
344 :
345 : !generate spherical harmonics
346 49296 : CALL harmonicsr(y, MATMUL(kghlp, cell%bmat), l)
347 :
348 3792 : lm = l**2
349 15168 : cdum = EXP(2*pi_const*img*dot_PRODUCT(kghlp, atoms%taual(:, ratom)))
350 3792 : IF (invsfct == 1) THEN
351 17072 : DO m = 1, 2*l + 1
352 13280 : lm = lm + 1
353 17072 : c_rot(m, igpt_lo, ilo, ratom, isym) = cdum*CONJG(y(lm))
354 : END DO
355 : ELSE
356 0 : DO m = 1, 2*l + 1
357 0 : lm = lm + 1
358 0 : c_rot(m, igpt_lo, ilo, ratom, isym) = cdum*CONJG(y(lm))
359 0 : c_rot(4*l + 3 - m, igpt_lo, ilo, ratom, isym) = (-1)**(m - 1)*(CONJG(cdum))*y(lm)
360 : END DO
361 : END IF
362 16608 : cfac(lapw%nv(jsp) + ilotot, isym) = EXP(-2*pi_const*img*(dot_PRODUCT(kghlp, trans(:, isym))))
363 : END DO
364 :
365 1440 : idum = invsfct*(2*l + 1)
366 :
367 4320 : allocate(ipiv(idum))
368 :
369 18512 : chelp(:idum, :idum) = c_lo(:idum, :idum, ilo, ratom)
370 :
371 33248 : CALL ZGESV(idum, idum, chelp(:idum, :idum), idum, ipiv, c_rot(:idum, :idum, ilo, ratom, isym), idum, ok)
372 :
373 1440 : IF (ok /= 0) call judft_error('symmetrizeh_new: failure zgesv')
374 2464 : deallocate(ipiv, y)
375 : END DO
376 : END IF
377 : END DO
378 : END DO
379 : END DO
380 :
381 : ! symmetrize local-orbital-apw-part of the matrix
382 24 : i = 0
383 24 : iatom = 0
384 60 : DO itype = 1, atoms%ntype
385 96 : DO ieq = 1, atoms%neq(itype)
386 36 : iatom = iatom + 1
387 72 : IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN
388 36 : IF (sym%invsat(iatom) == 0) invsfct = 1
389 36 : IF (sym%invsat(iatom) == 1) invsfct = 2
390 :
391 84 : DO ilo = 1, atoms%nlo(itype)
392 48 : l = atoms%llo(ilo, itype)
393 216 : DO igpt = 1, invsfct*(2*l + 1)
394 132 : i = i + 1
395 20200 : DO j = 1, lapw%nv(jsp)
396 : cdum = 0
397 : ic = 0
398 571444 : DO isym = 1, nsymop
399 : ic = ic + 1
400 551424 : iop = psym(isym)
401 551424 : ratom = map(isym, iatom)
402 551424 : IF (invsfct == 2) THEN
403 0 : IF (sym%invsat(ratom) == 2) THEN
404 0 : ratom = sym%invsatnr(ratom)
405 : END IF
406 : END IF
407 :
408 551424 : igpt_lo1 = lo_indx(ilo, ratom)
409 551424 : igpt_lo2 = igpt_lo1 + invsfct*2*l
410 551424 : IF (invsfct == 2) igpt_lo2 = igpt_lo2 + 1
411 551424 : igpt1 = pointer_apw(j, isym)
412 :
413 551424 : cdum2 = 0
414 551424 : ic1 = 0
415 2575584 : DO igpt2 = igpt_lo1, igpt_lo2
416 2024160 : ic1 = ic1 + 1
417 2575584 : IF (hmatTemp%l_real) THEN
418 732960 : cdum2 = cdum2 + CONJG(c_rot(ic1, igpt, ilo, ratom, isym))*hmatTemp%data_r(igpt1, lapw%nv(jsp) + igpt2)
419 : ELSE
420 1291200 : cdum2 = cdum2 + CONJG(c_rot(ic1, igpt, ilo, ratom, isym))*hmatTemp%data_c(igpt1, lapw%nv(jsp) + igpt2)
421 : END IF
422 : END DO
423 :
424 571444 : IF (iop <= sym%nop) THEN
425 422304 : cdum = cdum + cdum2*CONJG(cfac(lapw%nv(jsp) + i, isym))*cfac(j, isym)
426 : ELSE
427 129120 : cdum = cdum + CONJG(cdum2*CONJG(cfac(lapw%nv(jsp) + i, isym))*cfac(j, isym))
428 : END IF
429 : END DO
430 20152 : IF (hmat%l_real) THEN
431 9800 : hmat%data_r(j, lapw%nv(jsp) + i) = real(cdum)!/ic
432 9800 : hmat%data_r(lapw%nv(jsp) + i,j) = real(cdum)!/ic
433 : ELSE
434 10220 : hmat%data_c(j, lapw%nv(jsp) + i) = cdum!/ic
435 10220 : hmat%data_c(lapw%nv(jsp) + i,j) = CONJG(cdum)!/ic
436 : END IF
437 : END DO
438 : END DO
439 : END DO
440 : END IF
441 : END DO
442 : END DO
443 :
444 : !lo's - lo's
445 24 : i = 0
446 24 : iatom = 0
447 60 : DO itype = 1, atoms%ntype
448 96 : DO ieq = 1, atoms%neq(itype)
449 36 : iatom = iatom + 1
450 72 : IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN
451 36 : IF (sym%invsat(iatom) == 0) invsfct = 1
452 36 : IF (sym%invsat(iatom) == 1) invsfct = 2
453 :
454 84 : DO ilo = 1, atoms%nlo(itype)
455 48 : l = atoms%llo(ilo, itype)
456 216 : DO igpt = 1, invsfct*(2*l + 1)
457 132 : i = i + 1
458 132 : j = 0
459 132 : iatom1 = 0
460 396 : DO itype1 = 1, atoms%ntype
461 564 : DO ieq1 = 1, atoms%neq(itype1)
462 216 : iatom1 = iatom1 + 1
463 432 : IF ((sym%invsat(iatom1) == 0) .OR. (sym%invsat(iatom1) == 1)) THEN
464 216 : IF (sym%invsat(iatom1) == 0) invsfct1 = 1
465 216 : IF (sym%invsat(iatom1) == 1) invsfct1 = 2
466 :
467 480 : DO ilo1 = 1, atoms%nlo(itype1)
468 264 : l1 = atoms%llo(ilo1, itype1)
469 1368 : DO igpt1 = 1, invsfct1*(2*l1 + 1)
470 888 : j = j + 1
471 888 : IF (j > i) CYCLE
472 : cdum = 0
473 : ic = 0
474 14550 : DO isym = 1, nsymop
475 14040 : iop = psym(isym)
476 14040 : ratom = map(isym, iatom)
477 14040 : ratom1 = map(isym, iatom1)
478 :
479 14040 : IF (invsfct == 2) THEN
480 0 : IF (sym%invsat(ratom) == 2) THEN
481 0 : ratom = sym%invsatnr(ratom)
482 : END IF
483 : END IF
484 14040 : IF (invsfct1 == 2) THEN
485 0 : IF (sym%invsat(ratom1) == 2) THEN
486 0 : ratom1 = sym%invsatnr(ratom1)
487 : END IF
488 : END IF
489 :
490 14040 : igpt_lo1 = lo_indx(ilo, ratom)
491 14040 : igpt_lo2 = igpt_lo1 + invsfct*2*l
492 14040 : IF (invsfct == 2) igpt_lo2 = igpt_lo2 + 1
493 :
494 14040 : igpt1_lo1 = lo_indx(ilo1, ratom1)
495 14040 : igpt1_lo2 = igpt1_lo1 + invsfct1*2*l1
496 14040 : IF (invsfct1 == 2) igpt1_lo2 = igpt1_lo2 + 1
497 :
498 14040 : cdum2 = 0
499 14040 : ic1 = 0
500 71744 : DO igpt2 = igpt_lo1, igpt_lo2
501 57704 : ic1 = ic1 + 1
502 57704 : ic2 = 0
503 317096 : DO igpt3 = igpt1_lo1, igpt1_lo2
504 245352 : ic2 = ic2 + 1
505 303056 : IF (hmatTemp%l_real) THEN
506 : cdum2 = cdum2 + CONJG(c_rot(ic1, igpt, ilo, ratom, isym))* &
507 : hmatTemp%data_r(lapw%nv(jsp) + igpt3, lapw%nv(jsp) + igpt2)* &
508 36352 : c_rot(ic2, igpt1, ilo1, ratom1, isym)
509 : ELSE
510 : cdum2 = cdum2 + CONJG(c_rot(ic1, igpt, ilo, ratom, isym))* &
511 : hmatTemp%data_c(lapw%nv(jsp) + igpt3, lapw%nv(jsp) + igpt2)* &
512 209000 : c_rot(ic2, igpt1, ilo1, ratom1, isym)
513 : END IF
514 : END DO
515 : END DO
516 : ic = ic + 1
517 14550 : IF (iop <= sym%nop) THEN
518 9860 : cdum = cdum + cdum2*CONJG(cfac(lapw%nv(jsp) + i, isym))*cfac(lapw%nv(jsp) + j, isym)
519 : ELSE
520 4180 : cdum = cdum + CONJG(cdum2*CONJG(cfac(lapw%nv(jsp) + i, isym))*cfac(lapw%nv(jsp) + j, isym))
521 : END IF
522 : END DO
523 774 : IF (hmat%l_real) THEN
524 180 : hmat%data_r(lapw%nv(jsp) + j, lapw%nv(jsp) + i) = real(cdum)!/ic
525 : ELSE
526 330 : hmat%data_c(lapw%nv(jsp) + j, lapw%nv(jsp) + i) = cdum!/ic
527 : END IF
528 : END DO ! igpt_lo1
529 : END DO ! ilo1
530 : END IF
531 : END DO !ieq1
532 : END DO !itype1
533 : END DO ! igpt_lo
534 : END DO ! ilo
535 : END IF
536 : END DO !ieq
537 : END DO !itype
538 :
539 : END IF ! ANY(atoms%nlo.NE.0)
540 :
541 : CONTAINS
542 :
543 : ! Returns the spherical harmonics Y_lm(^rvec) for l = 0,...,ll in Y(1,...,(ll+1)**2).
544 3924 : SUBROUTINE harmonicsr(Y, rvec, ll)
545 : use m_judft
546 : use m_constants, only: CMPLX_NOT_INITALIZED
547 : IMPLICIT NONE
548 : INTEGER, INTENT(IN) :: ll
549 : REAL, INTENT(IN) :: rvec(:)
550 : COMPLEX, INTENT(INOUT) :: Y((ll + 1)**2)
551 : REAL :: stheta, ctheta, sphi, cphi, r, rvec1(3)
552 : INTEGER :: l, lm
553 : COMPLEX :: c
554 : COMPLEX, PARAMETER :: img = (0.0, 1.0)
555 :
556 25762 : Y = CMPLX_NOT_INITALIZED
557 3924 : Y(1) = 0.282094791773878
558 3924 : IF (ll == 0) RETURN
559 :
560 13352 : stheta = 0
561 13352 : ctheta = 0
562 13352 : sphi = 0
563 13352 : cphi = 0
564 13352 : r = norm2(rvec)
565 3338 : IF (r > 1e-16) THEN
566 13352 : rvec1 = rvec/r
567 3338 : ctheta = rvec1(3)
568 3338 : stheta = SQRT(rvec1(1)**2 + rvec1(2)**2)
569 3338 : IF (stheta > 1e-16) THEN
570 3014 : cphi = rvec1(1)/stheta
571 3014 : sphi = rvec1(2)/stheta
572 : END IF
573 : ELSE
574 0 : Y(2:) = 0.0
575 : RETURN
576 : END IF
577 :
578 : ! define Y,l,-l and Y,l,l
579 3338 : r = real(Y(1))
580 3338 : c = 1
581 8256 : DO l = 1, ll
582 4918 : r = r*stheta*SQRT(1.0 + 1.0/(2*l))
583 4918 : c = c*(cphi + img*sphi)
584 4918 : Y(l**2 + 1) = r*CONJG(c) ! l,-l
585 8256 : Y((l + 1)**2) = r*c*(-1)**l ! l,l
586 : END DO
587 :
588 : ! define Y,l,-l+1 and Y,l,l-1
589 3338 : Y(3) = 0.48860251190292*ctheta
590 4918 : DO l = 2, ll
591 1580 : r = SQRT(2.0*l + 1)*ctheta
592 1580 : Y(l**2 + 2) = r*Y((l - 1)**2 + 1) ! l,-l+1
593 4918 : Y(l*(l + 2)) = r*Y(l**2) ! l,l-1
594 : END DO
595 :
596 : ! define Y,l,m, |m|<l-1
597 4918 : DO l = 2, ll
598 1580 : lm = l**2 + 2
599 6498 : DO m = -l + 2, l - 2
600 1580 : lm = lm + 1
601 3160 : Y(lm) = SQRT((2.0*l + 1)/(l + m)/(l - m))*(SQRT(2.0*l - 1)*ctheta*Y(lm - 2*l) - SQRT((l + m - 1.0)*(l - m - 1)/(2*l - 3))*Y(lm - 4*l + 2))
602 : END DO
603 : END DO
604 : END SUBROUTINE harmonicsr
605 :
606 : END SUBROUTINE symmetrizeh
607 :
608 : END MODULE m_symmetrizeh
|
