Line data Source code
1 : MODULE m_types_hybdat
2 : use m_types_usdus
3 : use m_types_mat
4 : use m_types_coul
5 : use m_types_eigvec
6 : use m_io_matrix
7 : use m_judft
8 : use m_types_misc
9 : use m_types_fleurinput
10 : use m_types_mpi
11 : use m_constants
12 :
13 : #ifdef CPP_MPI
14 : use mpi
15 : #endif
16 : IMPLICIT NONE
17 : private
18 : TYPE,public:: t_hybdat
19 : COMPLEX, ALLOCATABLE :: stepfunc(:, :, :)
20 : INTEGER, ALLOCATABLE :: lmaxc(:)
21 : INTEGER, ALLOCATABLE :: nbands(:,:) ! nkptf, jsp
22 : INTEGER, ALLOCATABLE :: nbasm(:)
23 : INTEGER, ALLOCATABLE :: nindxc(:, :)
24 : INTEGER, ALLOCATABLE :: nindxp1(:, :)
25 : INTEGER, ALLOCATABLE :: nobd(:, :)
26 : INTEGER, ALLOCATABLE :: pntgpt(:, :, :, :)
27 : INTEGER, ALLOCATABLE :: pntgptd(:)
28 : INTEGER :: eig_id = -1
29 : INTEGER :: lmaxcd, maxindxc
30 : INTEGER :: maxfac
31 : INTEGER :: maxlmindx = -1
32 : INTEGER :: nbasp = -1
33 : integer :: max_q = -1
34 : LOGICAL :: l_addhf = .false.
35 : LOGICAL :: l_calhf = .false.
36 : LOGICAL :: l_print_iob_splitting = .True.
37 : LOGICAL :: l_subvxc = .false.
38 : REAL, ALLOCATABLE :: bas1(:, :, :, :), bas2(:, :, :, :)
39 : REAL, ALLOCATABLE :: bas1_MT(:, :, :), drbas1_MT(:, :, :)
40 : REAL, ALLOCATABLE :: core1(:, :, :, :), core2(:, :, :, :)
41 : REAL, ALLOCATABLE :: div_vv(:, :, :)
42 : REAL, ALLOCATABLE :: eig_c(:, :, :)
43 : REAL, ALLOCATABLE :: gauntarr(:, :, :, :, :, :)
44 : REAL, ALLOCATABLE :: gridf(:, :)
45 : REAL, ALLOCATABLE :: prodm(:, :, :, :)
46 : REAL, ALLOCATABLE :: sfac(:), fac(:)
47 : ! coulomb matrix stuff
48 : type(t_coul), allocatable :: coul(:)
49 :
50 : type(t_usdus) :: usdus
51 : class(t_mat), allocatable :: v_x(:, :) ! nkpt, jsp
52 : type(t_eigvec), allocatable :: zmat(:, :) ! nkpt, jsp
53 : type(t_results) :: results
54 : contains
55 : procedure :: set_stepfunction => set_stepfunction
56 : procedure :: free => free_hybdat
57 : procedure :: allocate => allocate_hybdat
58 : procedure :: set_nobd => set_nobd_hybdat
59 : procedure :: set_nbands => set_nbands_hybdat
60 : procedure :: set_maxlmindx => set_maxlmindx_hybdat
61 : END TYPE t_hybdat
62 : public:: gptnorm
63 : contains
64 0 : subroutine set_maxlmindx_hybdat(hybdat, atoms, num_radfun_per_l)
65 : implicit none
66 : class(t_hybdat), intent(inout) :: hybdat
67 : type(t_atoms), intent(in) :: atoms
68 : integer, intent(in) :: num_radfun_per_l(0:,:)
69 :
70 : integer :: itype, l, summe
71 :
72 0 : hybdat%maxlmindx = 0
73 0 : do itype = 1,atoms%ntype
74 0 : summe = 0
75 0 : do l = 0, atoms%lmax(itype)
76 0 : summe = summe + (2*l + 1) * num_radfun_per_l(l, itype)
77 : enddo
78 0 : hybdat%maxlmindx = max(hybdat%maxlmindx, summe)
79 : enddo
80 0 : end subroutine set_maxlmindx_hybdat
81 :
82 28 : subroutine set_nobd_hybdat(hybdat, fi, results)
83 : implicit none
84 : class(t_hybdat), intent(inout) :: hybdat
85 : type(t_fleurinput), intent(in) :: fi
86 : type(t_results), intent(in) :: results
87 :
88 : integer :: jsp, nk, i
89 :
90 28 : if (.not. allocated(hybdat%nobd)) then
91 24 : allocate (hybdat%nobd(fi%kpts%nkptf, fi%input%jspins))
92 : endif
93 388 : hybdat%nobd = 0
94 :
95 68 : do jsp = 1, fi%input%jspins
96 188 : DO nk = 1, fi%kpts%nkpt
97 8560 : DO i = 1, results%neig(nk, jsp)
98 8520 : IF (results%w_iks(i, nk, jsp) > 0.0) then
99 1116 : hybdat%nobd(nk, jsp) = hybdat%nobd(nk, jsp) + 1
100 : endif
101 : END DO
102 : enddo
103 : enddo
104 :
105 68 : do jsp = 1, fi%input%jspins
106 388 : DO nk = 1, fi%kpts%nkptf
107 320 : i = fi%kpts%bkp(nk)
108 360 : hybdat%nobd(nk, jsp) = hybdat%nobd(i, jsp)
109 : END DO
110 : enddo
111 28 : end subroutine set_nobd_hybdat
112 :
113 28 : subroutine set_nbands_hybdat(hybdat, fi, fmpi, results)
114 : implicit none
115 : class(t_hybdat), intent(inout) :: hybdat
116 : type(t_fleurinput), intent(in) :: fi
117 : type(t_mpi), intent(in) :: fmpi
118 : type(t_results), intent(in) :: results
119 :
120 :
121 : integer :: i, j, nk, jsp
122 :
123 28 : INTEGER :: degenerat(merge(fi%input%neig*2,fi%input%neig,fi%noco%l_soc) + 1, fi%kpts%nkpt)
124 : !determine degenerate states at each k-point
125 : !
126 : ! degenerat(i) =1 band i is not degenerat ,
127 : ! degenerat(i) =j band i has j-1 degenart states ( i, i+1, ..., i+j)
128 : ! degenerat(i) =0 band i is degenerat, but is not the lowest band
129 : ! of the group of degenerate states
130 :
131 28 : call timestart("degenerate treatment")
132 28 : IF (fmpi%irank == 0) THEN
133 14 : WRITE (oUnit, *)
134 14 : WRITE (oUnit, '(A)') " k-point | number of occupied bands | maximal number of bands"
135 : END IF
136 :
137 68 : do jsp = 1,fi%input%jspins
138 8680 : degenerat = 1
139 160 : DO nk = 1, fi%kpts%nkpt
140 8520 : DO i = 1, results%neig(nk, jsp)
141 307920 : DO j = i + 1, results%neig(nk, jsp)
142 307800 : IF (ABS(results%eig(i, nk, jsp) - results%eig(j, nk, jsp)) < 1E-07) THEN !0.015
143 4248 : degenerat(i, nk) = degenerat(i, nk) + 1
144 : END IF
145 : END DO
146 : END DO
147 :
148 8520 : DO i = 1, results%neig(nk, jsp)
149 11648 : IF ((degenerat(i, nk) /= 1) .OR. (degenerat(i, nk) /= 0)) degenerat(i + 1:i + degenerat(i, nk) - 1, nk) = 0
150 : END DO
151 :
152 : ! set the size of the exchange matrix in the space of the wavefunctions
153 :
154 120 : hybdat%nbands(nk,jsp) = fi%hybinp%bands1
155 120 : IF (hybdat%nbands(nk,jsp) > results%neig(nk, jsp)) THEN
156 0 : IF (fmpi%irank == 0) THEN
157 0 : WRITE (*, *) ' maximum for hybdat%nbands is', results%neig(nk, jsp)
158 0 : WRITE (*, *) ' increase energy window to obtain enough eigenvalues'
159 0 : WRITE (*, *) ' set hybdat%nbands equal to results%neig'
160 : END IF
161 0 : hybdat%nbands(nk,jsp) = results%neig(nk, jsp)
162 : END IF
163 :
164 220 : DO i = hybdat%nbands(nk,jsp) - 1, 1, -1
165 220 : IF ((degenerat(i, nk) >= 1) .AND. (degenerat(i, nk) + i - 1 /= hybdat%nbands(nk,jsp))) THEN
166 120 : hybdat%nbands(nk,jsp) = i + degenerat(i, nk) - 1
167 120 : EXIT
168 : END IF
169 : END DO
170 :
171 160 : IF (hybdat%nobd(nk, jsp) > hybdat%nbands(nk,jsp)) THEN
172 0 : WRITE (*, *) 'k-point: ', nk
173 0 : WRITE (*, *) 'number of bands: ', hybdat%nbands(nk,jsp)
174 0 : WRITE (*, *) 'number of occupied bands: ', hybdat%nobd(nk, jsp)
175 0 : CALL judft_warn("More occupied bands than total no of bands!?")
176 0 : hybdat%nbands(nk,jsp) = hybdat%nobd(nk, jsp)
177 : END IF
178 : END DO
179 :
180 : ! spread nbands from IBZ to whole BZ
181 268 : DO nk = fi%kpts%nkpt+1, fi%kpts%nkptf
182 200 : i = fi%kpts%bkp(nk)
183 240 : hybdat%nbands(nk,jsp) = hybdat%nbands(i,jsp)
184 : END DO
185 : enddo
186 28 : call timestop("degenerate treatment")
187 28 : end subroutine set_nbands_hybdat
188 :
189 12 : subroutine allocate_hybdat(hybdat, fi, num_radfun_per_l)
190 :
191 : implicit none
192 : class(t_hybdat), intent(inout) :: hybdat
193 : type(t_fleurinput), intent(in) :: fi
194 : integer, intent(in) :: num_radfun_per_l(:, :)
195 : integer :: ok(12)
196 :
197 156 : ok = -1
198 : allocate (hybdat%lmaxc(fi%atoms%ntype), &
199 56 : stat=ok(1), source=0)
200 : allocate (hybdat%bas1(fi%atoms%jmtd, maxval(num_radfun_per_l), 0:fi%atoms%lmaxd, fi%atoms%ntype), &
201 486024 : stat=ok(2), source=0.0)
202 : allocate (hybdat%bas2(fi%atoms%jmtd, maxval(num_radfun_per_l), 0:fi%atoms%lmaxd, fi%atoms%ntype), &
203 486024 : stat=ok(3), source=0.0)
204 : allocate (hybdat%bas1_MT(maxval(num_radfun_per_l), 0:fi%atoms%lmaxd, fi%atoms%ntype), &
205 1080 : stat=ok(4), source=0.0)
206 : allocate (hybdat%drbas1_MT(maxval(num_radfun_per_l), 0:fi%atoms%lmaxd, fi%atoms%ntype), &
207 1080 : stat=ok(5), source=0.0)
208 :
209 : ! core allocs
210 : allocate (hybdat%nindxc(0:hybdat%lmaxcd, fi%atoms%ntype), &
211 108 : stat=ok(6), source=0)
212 : allocate (hybdat%core1(fi%atoms%jmtd, hybdat%maxindxc, 0:hybdat%lmaxcd, fi%atoms%ntype), &
213 77348 : stat=ok(7), source=0.0)
214 : allocate (hybdat%core2(fi%atoms%jmtd, hybdat%maxindxc, 0:hybdat%lmaxcd, fi%atoms%ntype), &
215 77336 : stat=ok(8), source=0.0)
216 : allocate (hybdat%eig_c(hybdat%maxindxc, 0:hybdat%lmaxcd, fi%atoms%ntype), &
217 216 : stat=ok(9), source=0.0)
218 :
219 316 : allocate (hybdat%fac(0:hybdat%maxfac), stat=ok(10), source=0.0)
220 304 : allocate (hybdat%sfac(0:hybdat%maxfac), stat=ok(11), source=0.0)
221 :
222 : ALLOCATE (hybdat%gauntarr(2, 0:fi%atoms%lmaxd, 0:fi%atoms%lmaxd, 0:maxval(fi%hybinp%lcutm1), &
223 : -fi%atoms%lmaxd:fi%atoms%lmaxd, -maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1)), &
224 2968072 : stat=ok(12), source=0.0)
225 :
226 156 : if (any(ok /= 0)) then
227 0 : write (*, *) "allocation of hybdat failed. Error in array no.:"
228 0 : write (*, *) maxloc(abs(ok))
229 0 : call juDFT_error("allocation of hybdat failed. Error in array no is the outfile")
230 : endif
231 12 : end subroutine allocate_hybdat
232 :
233 12 : subroutine free_hybdat(hybdat)
234 : implicit none
235 : class(t_hybdat), intent(inout) :: hybdat
236 :
237 12 : if (allocated(hybdat%lmaxc)) deallocate (hybdat%lmaxc)
238 12 : if (allocated(hybdat%bas1)) deallocate (hybdat%bas1)
239 12 : if (allocated(hybdat%bas2)) deallocate (hybdat%bas2)
240 12 : if (allocated(hybdat%bas1_MT)) deallocate (hybdat%bas1_MT)
241 12 : if (allocated(hybdat%drbas1_MT)) deallocate (hybdat%drbas1_MT)
242 12 : if (allocated(hybdat%nindxc)) deallocate (hybdat%nindxc)
243 12 : if (allocated(hybdat%core1)) deallocate (hybdat%core1)
244 12 : if (allocated(hybdat%core2)) deallocate (hybdat%core2)
245 12 : if (allocated(hybdat%eig_c)) deallocate (hybdat%eig_c)
246 12 : if (allocated(hybdat%fac)) deallocate (hybdat%fac)
247 12 : if (allocated(hybdat%sfac)) deallocate (hybdat%sfac)
248 12 : if (allocated(hybdat%gauntarr)) deallocate (hybdat%gauntarr)
249 12 : end subroutine free_hybdat
250 :
251 0 : subroutine set_stepfunction(hybdat, cell, atoms, g, svol)
252 :
253 : implicit none
254 : class(t_hybdat), INTENT(INOUT) :: hybdat
255 : type(t_cell), INTENT(in) :: cell
256 : type(t_atoms), INTENT(in) :: atoms
257 : integer, INTENT(in) :: g(3)
258 : real, INTENT(in) :: svol
259 : integer :: i, j, k, ok
260 :
261 0 : if (.not. allocated(hybdat%stepfunc)) then
262 0 : call timestart("setup stepfunction")
263 0 : ALLOCATE (hybdat%stepfunc(-g(1):g(1), -g(2):g(2), -g(3):g(3)), stat=ok)
264 0 : IF (ok /= 0) then
265 0 : call juDFT_error('wavefproducts_inv5: error allocation stepfunc')
266 : endif
267 :
268 0 : DO i = -g(1), g(1)
269 0 : DO j = -g(2), g(2)
270 0 : DO k = -g(3), g(3)
271 0 : hybdat%stepfunc(i, j, k) = stepfunction(cell, atoms, [i, j, k])/svol
272 : END DO
273 : END DO
274 : END DO
275 0 : call timestop("setup stepfunction")
276 : endif
277 :
278 0 : end subroutine set_stepfunction
279 :
280 : !private subroutine
281 0 : FUNCTION stepfunction(cell, atoms, g)
282 :
283 : IMPLICIT NONE
284 :
285 : TYPE(t_cell), INTENT(IN) :: cell
286 : TYPE(t_atoms), INTENT(IN) :: atoms
287 :
288 : INTEGER, INTENT(IN) :: g(3)
289 : COMPLEX :: stepfunction !Is real in inversion case
290 : REAL :: gnorm, gnorm3, r, fgr
291 : INTEGER :: itype, ieq, icent
292 :
293 0 : gnorm = gptnorm(g, cell%bmat)
294 0 : gnorm3 = gnorm**3
295 0 : IF (abs(gnorm) < 1e-12) THEN
296 0 : stepfunction = 1
297 0 : DO itype = 1, atoms%ntype
298 0 : stepfunction = stepfunction - atoms%neq(itype)*atoms%volmts(itype)/cell%omtil
299 : END DO
300 : ELSE
301 0 : stepfunction = 0
302 0 : icent = 0
303 0 : DO itype = 1, atoms%ntype
304 0 : r = gnorm*atoms%rmt(itype)
305 0 : fgr = fpi_const*(sin(r) - r*cos(r))/gnorm3/cell%omtil
306 0 : DO ieq = 1, atoms%neq(itype)
307 0 : icent = icent + 1
308 0 : stepfunction = stepfunction - fgr*exp(-cmplx(0., tpi_const*dot_product(atoms%taual(:, icent), g)))
309 : ENDDO
310 : ENDDO
311 : ENDIF
312 :
313 0 : END FUNCTION stepfunction
314 :
315 988068 : PURE FUNCTION gptnorm(gpt, bmat)
316 : IMPLICIT NONE
317 : REAL :: gptnorm
318 : INTEGER, INTENT(IN) :: gpt(3)
319 : REAL, INTENT(IN) :: bmat(3, 3)
320 :
321 18773292 : gptnorm = norm2(matmul(gpt(:), bmat(:, :)))
322 :
323 988068 : END FUNCTION gptnorm
324 :
325 : subroutine calc_matrix_slots(l_real, mtx_dim, slots_per_mtx, col_in_slot)
326 : implicit none
327 : logical, intent(in) :: l_real
328 : integer, intent(in) :: mtx_dim
329 : integer, intent(out) :: slots_per_mtx, col_in_slot
330 :
331 : integer(8) :: mtx_size, type_size
332 : integer(8), parameter :: max_bytes = huge(slots_per_mtx) - 1
333 : integer :: i
334 : type_size = merge(8, 16, l_real)
335 :
336 : ! avoid int32 overflow
337 : mtx_size = type_size*mtx_dim
338 : mtx_size = mtx_size*mtx_dim
339 :
340 : i = 1
341 : do while ((mtx_size/i >= max_bytes) .or. mod(mtx_dim, i) /= 0)
342 : i = i + 1
343 : enddo
344 :
345 : slots_per_mtx = i
346 : col_in_slot = mtx_dim/i
347 : end subroutine calc_matrix_slots
348 0 : END MODULE m_types_hybdat
|
