Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : ! This module generates the little group of k and the extended irr. !
3 : ! BZ. Furthermore it calculates the irr. representation !
4 : ! !
5 : ! P(R,T)\phi_n,k = \sum_{n'} rep_v(n',n) *\phi_n',k !
6 : ! where !
7 : ! P is an element of the little group of k !
8 : ! n' runs over group of degenerat states belonging to n. !
9 : ! !
10 : ! M.Betzinger (09/07) !
11 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
12 : MODULE m_symm_hf
13 :
14 : use m_judft
15 : USE m_types
16 : USE m_types_hybdat
17 : USE m_constants
18 : USE m_util
19 : USE m_intgrf
20 : USE m_io_hybrid
21 : #ifdef CPP_MPI
22 : use mpi
23 : #endif
24 : CONTAINS
25 :
26 24 : SUBROUTINE symm_hf_init(fi, nk, nsymop, rrot, psym)
27 : use m_juDFT
28 : IMPLICIT NONE
29 :
30 : type(t_fleurinput), intent(in) :: fi
31 : INTEGER, INTENT(IN) :: nk
32 : INTEGER, INTENT(INOUT) :: nsymop
33 : INTEGER, INTENT(INOUT) :: rrot(:, :, :) ! 3,3,fi%sym%nsym
34 : INTEGER, INTENT(INOUT) :: psym(:) ! Note: psym is only filled up to index nsymop
35 :
36 : INTEGER :: i
37 : REAL :: rotkpt(3)
38 :
39 24 : CALL timestart("symm_hf_init")
40 24 : nsymop = 0
41 : ! calculate rotations in reciprocal space
42 1176 : DO i = 1, fi%sym%nsym
43 1176 : IF(i <= fi%sym%nop) THEN
44 13104 : rrot(:, :, i) = transpose(fi%sym%mrot(:, :, fi%sym%invtab(i)))
45 : ELSE
46 1872 : rrot(:, :, i) = -rrot(:, :, i - fi%sym%nop)
47 : END IF
48 : END DO
49 :
50 : ! determine little group of k., i.e. those symmetry operations
51 : ! which keep bk(:,nk) invariant
52 : ! nsymop :: number of such symmetry-operations
53 : ! psym :: points to the symmetry-operation
54 :
55 1176 : psym = 0
56 : nsymop = 0
57 1176 : DO i = 1, fi%sym%nsym
58 28800 : rotkpt = matmul(rrot(:, :, i), fi%kpts%bkf(:, nk))
59 :
60 : !transfer rotkpt into BZ
61 4608 : rotkpt = fi%kpts%to_first_bz(rotkpt)
62 :
63 : !check if rotkpt is identical to bk(:,nk)
64 4632 : IF(maxval(abs(rotkpt - fi%kpts%bkf(:, nk))) <= 1E-07) THEN
65 720 : nsymop = nsymop + 1
66 720 : psym(nsymop) = i
67 : END IF
68 : END DO
69 24 : CALL timestop("symm_hf_init")
70 24 : END SUBROUTINE symm_hf_init
71 :
72 72 : SUBROUTINE symm_hf(fi, nk, hybdat, results, submpi, eig_irr, mpdata, cmt, &
73 24 : rrot, nsymop, psym, n_q, parent, nsest, indx_sest, jsp)
74 :
75 : USE m_olap
76 : USE m_trafo
77 : use m_calc_cmt
78 : use m_juDFT
79 :
80 : IMPLICIT NONE
81 :
82 : type(t_fleurinput), intent(in) :: fi
83 : TYPE(t_hybdat), INTENT(IN) :: hybdat
84 : type(t_results),intent(in) :: results
85 : type(t_hybmpi), intent(in) :: submpi
86 : TYPE(t_mpdata), intent(in) :: mpdata
87 :
88 : ! - scalars -
89 : INTEGER, INTENT(IN) :: nk, jsp
90 : INTEGER, INTENT(IN) :: nsymop
91 :
92 : ! - arrays -
93 : COMPLEX , intent(in) :: cmt(hybdat%nbands(nk,jsp), hybdat%maxlmindx, fi%atoms%nat)
94 : INTEGER, INTENT(IN) :: rrot(:, :, :)
95 : INTEGER, INTENT(IN) :: psym(:)
96 : INTEGER, INTENT(INOUT) :: parent(fi%kpts%nkptf)
97 : INTEGER, INTENT(INOUT) :: nsest(hybdat%nbands(nk,jsp))
98 : INTEGER, INTENT(INOUT) :: indx_sest(hybdat%nbands(nk,jsp), hybdat%nbands(nk,jsp))
99 : INTEGER, ALLOCATABLE, INTENT(INOUT) :: n_q(:)
100 :
101 : REAL, INTENT(IN) :: eig_irr(:, :)
102 :
103 : ! - local scalars -
104 : INTEGER :: ikpt, iop, isym, iisym, m
105 : INTEGER :: itype, ieq, iatom, ratom, ierr
106 : INTEGER :: iband1, iband2, iatom0
107 : INTEGER :: i, j, ic, ic1, ic2
108 : INTEGER :: ok, ld_olapmt, ld_cmthlp, ld_tmp, ld_wavefolap
109 : INTEGER :: l, lm
110 : INTEGER :: n1, n2, nn
111 : INTEGER :: ndb1, ndb2
112 : INTEGER :: nrkpt
113 : INTEGER :: maxndb, nddb
114 :
115 : REAL :: tolerance, pi
116 :
117 : COMPLEX :: cdum
118 : COMPLEX, PARAMETER :: img = (0.0, 1.0)
119 :
120 : ! - local arrays -
121 24 : INTEGER :: neqvkpt(fi%kpts%nkptf)
122 24 : INTEGER :: list(fi%kpts%nkptf)
123 24 : INTEGER :: degenerat(results%neig(nk, jsp))
124 :
125 : REAL :: rotkpt(3), g(3)
126 24 : complex, ALLOCATABLE :: olapmt(:, :, :, :)
127 :
128 24 : COMPLEX, ALLOCATABLE :: carr(:), wavefolap(:, :), tmp(:,:), carr_tmp(:)
129 24 : COMPLEX, ALLOCATABLE :: cmthlp(:, :)
130 24 : LOGICAL, ALLOCATABLE :: symequivalent(:, :)
131 :
132 24 : CALL timestart("symm_hf")
133 67096 : parent = 0; nsest = 0; indx_sest = 0;
134 :
135 : ! determine extented irreducible BZ of k ( EIBZ(k) ), i.e.
136 : ! those k-points, which can generate the whole BZ by
137 : ! applying the symmetry operations of the little group of k
138 24 : call timestart("calc EIBZ")
139 216 : neqvkpt = 0
140 :
141 216 : DO i = 1, fi%kpts%nkptf
142 216 : list(i) = i - 1
143 : END DO
144 :
145 192 : DO ikpt = 2, fi%kpts%nkptf
146 1292 : DO iop = 1, nsymop
147 :
148 30300 : rotkpt = matmul(rrot(:, :, psym(iop)), fi%kpts%bkf(:, ikpt))
149 :
150 : !determine number of rotkpt
151 1212 : nrkpt = fi%kpts%get_nk(rotkpt)
152 1212 : IF(nrkpt == 0) call judft_error('symm: Difference vector not found !')
153 :
154 1212 : IF(list(nrkpt) /= 0) THEN
155 168 : list(nrkpt) = 0
156 168 : neqvkpt(ikpt) = neqvkpt(ikpt) + 1
157 168 : parent(nrkpt) = ikpt
158 : END IF
159 5480 : IF(all(list == 0)) EXIT
160 :
161 : END DO
162 : END DO
163 :
164 : ! for the Gamma-point holds:
165 24 : parent(1) = 1
166 24 : neqvkpt(1) = 1
167 24 : call timestop("calc EIBZ")
168 :
169 : ! determine the factor n_q, that means the number of symmetrie operations of the little group of bk(:,nk)
170 : ! which keep q (in EIBZ) invariant
171 24 : call timestart("calc n_q")
172 24 : IF(ALLOCATED(n_q)) DEALLOCATE(n_q)
173 160 : allocate(n_q(fi%kpts%EIBZ(nk)%nkpt), source=0)
174 :
175 112 : ic = 0
176 112 : n_q = 0
177 216 : DO ikpt = 1, fi%kpts%nkptf
178 216 : IF(parent(ikpt) == ikpt) THEN
179 88 : ic = ic + 1
180 2424 : DO iop = 1, nsymop
181 2336 : isym = psym(iop)
182 58400 : rotkpt = matmul(rrot(:, :, isym), fi%kpts%bkf(:, ikpt))
183 :
184 : !transfer rotkpt into BZ
185 9344 : rotkpt = fi%kpts%to_first_bz(rotkpt)
186 :
187 : !check if rotkpt is identical to bk(:,ikpt)
188 9432 : IF(maxval(abs(rotkpt - fi%kpts%bkf(:, ikpt))) <= 1E-06) THEN
189 1576 : n_q(ic) = n_q(ic) + 1
190 : END IF
191 : END DO
192 : END IF
193 : END DO
194 24 : IF(ic /= fi%kpts%EIBZ(nk)%nkpt) call judft_error('symm: failure EIBZ')
195 24 : call timestop("calc n_q")
196 :
197 : ! calculate degeneracy:
198 : ! degenerat(i) = 1 state i is not degenerat,
199 : ! degenerat(i) = j state i has j-1 degenerat states at {i+1,...,i+j-1}
200 : ! degenerat(i) = 0 state i is degenerat
201 24 : call timestart("calculate degeneracy")
202 24 : tolerance = 1E-07 !0.00001
203 :
204 1644 : degenerat = 1
205 :
206 1252 : DO i = 1, hybdat%nbands(nk,jsp)
207 32826 : DO j = i + 1, hybdat%nbands(nk,jsp)
208 32802 : IF(abs(eig_irr(i, nk) - eig_irr(j, nk)) <= tolerance) THEN
209 600 : degenerat(i) = degenerat(i) + 1
210 : END IF
211 : END DO
212 : END DO
213 :
214 1644 : DO i = 1, results%neig(nk, jsp)
215 24 : IF(degenerat(i) /= 1 .or. degenerat(i) /= 0) THEN
216 2072 : degenerat(i + 1:i + degenerat(i) - 1) = 0
217 : END IF
218 : END DO
219 :
220 : ! maximal number of degenerate bands -> maxndb
221 : maxndb = maxval(degenerat)
222 :
223 : ! number of different degenerate bands/states
224 1644 : nddb = count(degenerat >= 1)
225 24 : call timestop("calculate degeneracy")
226 :
227 24 : call timestart("calc olapmt")
228 : IF(allocated(olapmt)) deallocate(olapmt)
229 528 : allocate(olapmt(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:fi%atoms%lmaxd, fi%atoms%ntype), stat=ok)
230 24 : IF(ok /= 0) call judft_error('symm: failure allocation olapmt')
231 4584 : olapmt = 0
232 :
233 60 : DO itype = 1, fi%atoms%ntype
234 408 : DO l = 0, fi%atoms%lmax(itype)
235 348 : nn = mpdata%num_radfun_per_l(l, itype)
236 1128 : DO n2 = 1, nn
237 2724 : DO n1 = 1, nn
238 : olapmt(n1, n2, l, itype) = intgrf( &
239 : hybdat%bas1(:, n1, l, itype)*hybdat%bas1(:, n2, l, itype) &
240 : + hybdat%bas2(:, n1, l, itype)*hybdat%bas2(:, n2, l, itype), &
241 1324824 : fi%atoms, itype, hybdat%gridf)
242 : END DO
243 : END DO
244 : END DO
245 : END DO
246 24 : call timestop("calc olapmt")
247 :
248 552 : allocate(wavefolap(hybdat%nbands(nk,jsp), hybdat%nbands(nk,jsp)), carr(maxval(mpdata%num_radfun_per_l)), stat=ok)
249 24 : IF(ok /= 0) call judft_error('symm: failure allocation wfolap/maxindx')
250 65628 : wavefolap = 0
251 :
252 96 : allocate(cmthlp(size(cmt,2), size(cmt,1) ), stat=ierr)
253 24 : if(ierr /= 0) call judft_error("can't alloc cmthlp")
254 :
255 480 : allocate(tmp(maxval(mpdata%num_radfun_per_l), hybdat%nbands(nk,jsp)), stat=ierr)
256 24 : if(ierr /= 0 ) call judft_error("cant't alloc tmp")
257 :
258 24 : call timestart("calc wavefolap")
259 24 : ld_olapmt = size(olapmt,1)
260 24 : ld_cmthlp = size(cmthlp,1)
261 24 : ld_tmp = size(tmp, 1)
262 24 : ld_wavefolap = size(wavefolap,1)
263 :
264 60 : do iatom = 1+submpi%rank, fi%atoms%nat, submpi%size
265 36 : itype = fi%atoms%itype(iatom)
266 36 : call timestart("transp cmthlp")
267 331352 : cmthlp = transpose(cmt(:,:,iatom))
268 36 : call timestop("transp cmthlp")
269 36 : lm = 0
270 408 : DO l = 0, fi%atoms%lmax(itype)
271 3780 : DO M = -l, l
272 3396 : nn = mpdata%num_radfun_per_l(l, itype)
273 :
274 : !call zgemm(transa, transb, m, n, k, alpha, a, lda,
275 : call zgemm("N", "N", nn, hybdat%nbands(nk,jsp), nn, cmplx_1, olapmt(1,1,l,itype), ld_olapmt, &
276 : ! b, ldb, beta, c, ldc)
277 3396 : cmthlp(lm+1,1), ld_cmthlp, cmplx_0, tmp, ld_tmp)
278 :
279 : !call zgemm(transa, transb, m, n, k, alpha, a, lda,
280 : call zgemm("C", "N", hybdat%nbands(nk,jsp), hybdat%nbands(nk,jsp), nn, cmplx_1, cmthlp(lm+1, 1), ld_cmthlp, &
281 : ! b, ldb, beta, c, ldc)
282 3396 : tmp, ld_tmp, cmplx_1, wavefolap, ld_wavefolap)
283 3744 : lm = lm + nn
284 : END DO
285 : END DO
286 : END DO
287 :
288 24 : deallocate(cmthlp)
289 : #ifdef CPP_MPI
290 24 : call timestart("allreduce wavefolap")
291 72 : call MPI_ALLREDUCE(MPI_IN_PLACE, wavefolap, size(wavefolap), MPI_DOUBLE_COMPLEX, MPI_SUM, submpi%comm, ierr)
292 24 : call timestop("allreduce wavefolap")
293 : #endif
294 24 : call timestop("calc wavefolap")
295 :
296 24 : call timestart("calc symmequivalent")
297 :
298 63700 : allocate(symequivalent(nddb, nddb), stat=ok, source=.False.)
299 24 : IF(ok /= 0) call judft_error('symm: failure allocation symequivalent')
300 :
301 : !$OMP PARALLEL DO default(none) private(iband1, ndb1, ic1, iband2, ndb2, ic2) &
302 24 : !$OMP shared(submpi, hybdat, degenerat, wavefolap, symequivalent, nk, jsp)
303 : DO iband1 = submpi%rank + 1, hybdat%nbands(nk,jsp), submpi%size
304 : ndb1 = degenerat(iband1)
305 : IF(ndb1 /= 0) then
306 : ic1 = count(degenerat(:iband1) /= 0)
307 : DO iband2 = 1, hybdat%nbands(nk,jsp)
308 : ndb2 = degenerat(iband2)
309 : IF(ndb2 /= 0) then
310 : ic2 = count(degenerat(:iband2) /= 0)
311 : IF(any(abs(wavefolap(iband1:iband1 + ndb1 - 1, &
312 : iband2:iband2 + ndb2 - 1)) > 1E-9)) THEN
313 : symequivalent(ic2, ic1) = .true.
314 : END IF
315 : endif
316 : END DO
317 : endif
318 : END DO
319 : !$OMP end parallel do
320 : #ifdef CPP_MPI
321 24 : call timestart("allreduce symequivalent")
322 72 : call MPI_ALLREDUCE(MPI_IN_PLACE, symequivalent, size(symequivalent), MPI_LOGICAL, MPI_LOR, submpi%comm, ierr)
323 24 : call timestop("allreduce symequivalent")
324 : #endif
325 24 : call timestop("calc symmequivalent")
326 : !
327 : ! generate index field which contain the band combinations (n1,n2),
328 : ! which are non zero
329 : !
330 24 : call timestart("calc bandcombos")
331 24 : ic1 = 0
332 65628 : indx_sest = 0
333 1252 : nsest = 0
334 1252 : DO iband1 = 1, hybdat%nbands(nk,jsp)
335 1228 : ndb1 = degenerat(iband1)
336 1228 : IF(ndb1 >= 1) ic1 = ic1 + 1
337 1228 : i = 0
338 1228 : DO WHILE(degenerat(iband1 - i) == 0)
339 : i = i + 1
340 : END DO
341 1228 : ndb1 = degenerat(iband1 - i)
342 1228 : ic2 = 0
343 65628 : DO iband2 = 1, hybdat%nbands(nk,jsp)
344 64376 : ndb2 = degenerat(iband2)
345 64376 : IF(ndb2 >= 1) ic2 = ic2 + 1
346 64376 : i = 0
347 64376 : DO WHILE(degenerat(iband2 - i) == 0)
348 : i = i + 1
349 : END DO
350 64376 : ndb2 = degenerat(iband2 - i)
351 : ! only upper triangular part
352 65604 : IF(symequivalent(ic2, ic1) .and. iband2 <= iband1) THEN
353 : ! IF( ndb1 .ne. ndb2 ) call judft_error('symm_hf: failure symequivalent')
354 8008 : nsest(iband1) = nsest(iband1) + 1
355 8008 : indx_sest(nsest(iband1), iband1) = iband2
356 : END IF
357 : END DO
358 : END DO
359 24 : call timestop("calc bandcombos")
360 :
361 : !
362 : ! calculate representations for core states
363 : ! (note that for a core state, these are proportional to the Wigner D matrices)
364 : !
365 : ! Definition of the Wigner rotation matrices
366 : !
367 : ! -1 l
368 : ! P(R) Y (r) = Y (R r) = sum D (R) Y (r)
369 : ! lm lm m' m m' lm'
370 : !
371 :
372 : pi = pimach()
373 :
374 24 : call timestart("calc core repr")
375 24 : IF(hybdat%lmaxcd > fi%atoms%lmaxd) then
376 0 : call judft_error('symm_hf: The very impropable case that hybdat%lmaxcd > fi%atoms%lmaxd occurs')
377 : endif
378 24 : iatom = 0
379 24 : iatom0 = 0
380 60 : DO itype = 1, fi%atoms%ntype
381 72 : DO ieq = 1, fi%atoms%neq(itype)
382 36 : iatom = iatom + 1
383 1096 : DO iop = 1, nsymop
384 1024 : isym = psym(iop)
385 1024 : IF(isym <= fi%sym%nop) THEN
386 : iisym = isym
387 : ELSE
388 : iisym = isym - fi%sym%nop
389 : END IF
390 :
391 1024 : ratom = fi%hybinp%map(iatom, isym)
392 25600 : rotkpt = matmul(rrot(:, :, isym), fi%kpts%bkf(:, nk))
393 : g = nint(rotkpt - fi%kpts%bkf(:, nk))
394 :
395 : cdum = exp(-2*pi*img*dot_product(rotkpt, fi%sym%tau(:, iisym)))* &
396 36 : exp(2*pi*img*dot_product(g, fi%atoms%taual(:, ratom)))
397 : END DO
398 : END DO
399 60 : iatom0 = iatom0 + fi%atoms%neq(itype)
400 : END DO
401 24 : call timestop("calc core repr")
402 :
403 24 : CALL timestop("symm_hf")
404 24 : END SUBROUTINE symm_hf
405 5724 : END MODULE m_symm_hf
|
