Line data Source code
1 : module m_types_mpdata
2 : implicit none
3 :
4 : type t_mpdata
5 : integer, allocatable :: g(:, :) ! (3, num_gpts)
6 : integer, allocatable :: n_g(:) ! (ik)
7 : integer, allocatable :: gptm_ptr(:, :) ! (ig, ik)
8 : integer, allocatable :: num_radbasfn(:, :) !(l,itype)
9 : real, allocatable :: radbasfn_mt(:, :, :, :) !(jri,n,l,itype)
10 : INTEGER, ALLOCATABLE :: num_radfun_per_l(:, :) !(l,itype)
11 : INTEGER :: max_indx_p_1 = -1
12 :
13 : integer, allocatable :: l1(:, :, :) !(n, l, itype)
14 : integer, allocatable :: l2(:, :, :) !(n, l, itype)
15 : integer, allocatable :: n1(:, :, :) !(n, l, itype)
16 : integer, allocatable :: n2(:, :, :) !(n, l, itype)
17 : CONTAINS
18 : procedure :: num_gpts => mpdata_num_gpts
19 : procedure :: gen_gvec => mpdata_gen_gvec
20 : procedure :: check_orthonormality => mpdata_check_orthonormality
21 : procedure :: check_radbasfn => mpdata_check_radbasfn
22 : procedure :: calc_olap_radbasfn => mpdata_calc_olap_radbasfn
23 : procedure :: filter_radbasfn => mpdata_filter_radbasfn
24 : procedure :: trafo_to_orthonorm_bas => mpdata_trafo_to_orthonorm_bas
25 : procedure :: add_l0_fun => mpdata_add_l0_fun
26 : procedure :: reduce_linear_dep => mpdata_reduce_linear_dep
27 : procedure :: normalize => mpdata_normalize
28 : procedure :: set_nl => mpdata_set_nl
29 : procedure :: free => mpdata_free
30 : procedure :: init => mpdata_init
31 : procedure :: set_num_radfun_per_l => set_num_radfun_per_l_mpdata
32 : procedure :: set_max_indx_p_1 => set_max_indx_p_1
33 : !generic :: write(unformatted) => write_mpdata
34 : end type t_mpdata
35 : contains
36 30 : function mpdata_num_gpts(mpdata)
37 : implicit NONE
38 : class(t_mpdata), intent(in) :: mpdata
39 : integer :: mpdata_num_gpts
40 :
41 30 : mpdata_num_gpts = size(mpdata%g, dim=2)
42 30 : end function mpdata_num_gpts
43 :
44 12 : subroutine mpdata_gen_gvec(mpdata, mpinp, cell, kpts, mpi)
45 : use m_judft
46 : use m_types_setup
47 : use m_types_kpts
48 : use m_types_mpi
49 : use m_types_mpinp
50 : USE m_constants
51 : use m_intgrf, only: intgrf_init, intgrf
52 : use m_sort
53 : implicit NONE
54 : class(t_mpdata), intent(inout) :: mpdata
55 : type(t_mpinp), intent(in) :: mpinp
56 : type(t_cell), intent(in) :: cell
57 : type(t_kpts), intent(in) :: kpts
58 : type(t_mpi), intent(in) :: mpi
59 :
60 : integer :: i, n, n1, n2, divconq
61 : integer :: x, y, z, ikpt, igpt
62 : integer :: g(3)
63 : real :: longest_k, kvec(3)
64 : logical :: l_found_new_gpt, l_found_kg_in_sphere
65 :
66 12 : integer, allocatable :: unsrt_pgptm(:, :) ! unsorted pointers to g vectors
67 12 : real, allocatable :: length_kg(:, :) ! length of the vectors k + G
68 12 : integer, allocatable :: ptr(:)
69 :
70 36 : allocate(mpdata%n_g(kpts%nkptf))
71 :
72 108 : mpdata%n_g = 0
73 12 : i = 0
74 12 : n = -1
75 :
76 492 : longest_k = MAXVAL([(norm2(MATMUL(kpts%bkf(:, ikpt), cell%bmat)), ikpt=1, kpts%nkptf)])
77 :
78 : ! a first run for the determination of the dimensions of the fields g,pgptm
79 :
80 : do
81 128 : n = n + 1
82 128 : l_found_new_gpt = .FALSE.
83 1576 : do x = -n, n
84 1448 : n1 = n - ABS(x)
85 13000 : do y = -n1, n1
86 11424 : n2 = n1 - ABS(y)
87 33068 : do z = -n2, n2, MAX(2*n2, 1)
88 80784 : g = [x, y, z]
89 383724 : if((norm2(MATMUL(g, cell%bmat)) - longest_k) > mpinp%g_cutoff) CYCLE
90 3092 : l_found_kg_in_sphere = .FALSE.
91 39252 : do ikpt = 1, kpts%nkptf
92 497812 : if(norm2(MATMUL(kpts%bkf(:, ikpt) + g, cell%bmat)) <= mpinp%g_cutoff) THEN
93 13244 : if(.NOT. l_found_kg_in_sphere) THEN
94 2368 : i = i + 1
95 2368 : l_found_kg_in_sphere = .TRUE.
96 : END if
97 :
98 13244 : mpdata%n_g(ikpt) = mpdata%n_g(ikpt) + 1
99 13244 : l_found_new_gpt = .TRUE.
100 : END if
101 : enddo ! k-loop
102 : enddo
103 : enddo
104 : enddo
105 128 : if(.NOT. l_found_new_gpt) EXIT
106 : enddo
107 :
108 36 : allocate(mpdata%g(3, i)) ! i = gptmd
109 144 : allocate(mpdata%gptm_ptr(maxval(mpdata%n_g), kpts%nkptf))
110 :
111 : ! allocate and initialize arrays needed for G vector ordering
112 144 : allocate(unsrt_pgptm(maxval(mpdata%n_g), kpts%nkptf))
113 144 : allocate(length_kG(maxval(mpdata%n_g), kpts%nkptf))
114 :
115 9484 : mpdata%g = 0
116 14380 : mpdata%gptm_ptr = 0
117 108 : mpdata%n_g = 0
118 :
119 108 : i = 0
120 108 : n = -1
121 :
122 14380 : length_kG = 0
123 14380 : unsrt_pgptm = 0
124 :
125 : do
126 128 : n = n + 1
127 128 : l_found_new_gpt = .FALSE.
128 1576 : do x = -n, n
129 1448 : n1 = n - ABS(x)
130 13000 : do y = -n1, n1
131 11424 : n2 = n1 - ABS(y)
132 33068 : do z = -n2, n2, MAX(2*n2, 1)
133 80784 : g = [x, y, z]
134 383724 : if((norm2(MATMUL(g, cell%bmat)) - longest_k) > mpinp%g_cutoff) CYCLE
135 3092 : l_found_kg_in_sphere = .FALSE.
136 39252 : do ikpt = 1, kpts%nkptf
137 98944 : kvec = kpts%bkf(:, ikpt)
138 :
139 473076 : if(norm2(MATMUL(kvec + g, cell%bmat)) <= mpinp%g_cutoff) THEN
140 13244 : if(.NOT. l_found_kg_in_sphere) THEN
141 2368 : i = i + 1
142 9472 : mpdata%g(:, i) = g
143 : l_found_kg_in_sphere = .TRUE.
144 : END if
145 :
146 13244 : mpdata%n_g(ikpt) = mpdata%n_g(ikpt) + 1
147 13244 : l_found_new_gpt = .TRUE.
148 :
149 : ! Position of the vector is saved as pointer
150 13244 : unsrt_pgptm(mpdata%n_g(ikpt), ikpt) = i
151 : ! Save length of vector k + G for array sorting
152 251636 : length_kG(mpdata%n_g(ikpt), ikpt) = norm2(MATMUL(kvec + g, cell%bmat))
153 : END if
154 : enddo
155 : enddo
156 : enddo
157 : enddo
158 128 : if(.NOT. l_found_new_gpt) EXIT
159 : enddo
160 :
161 : ! Sort pointers in array, so that shortest |k+G| comes first
162 108 : do ikpt = 1, kpts%nkptf
163 288 : allocate(ptr(mpdata%n_g(ikpt)))
164 : ! create pointers which correspond to a sorted array (make algo stable)
165 26680 : call sort(ptr, length_kG(1:mpdata%n_g(ikpt), ikpt), [(1.0*i, i=1,mpdata%n_g(ikpt))])
166 : ! rearrange old pointers
167 13340 : do igpt = 1, mpdata%n_g(ikpt)
168 13340 : mpdata%gptm_ptr(igpt, ikpt) = unsrt_pgptm(ptr(igpt), ikpt)
169 : enddo
170 108 : deallocate(ptr)
171 : enddo
172 :
173 12 : if(mpi%irank == 0) THEN
174 6 : WRITE(oUnit, '(/A)') 'Mixed basis'
175 6 : WRITE(oUnit, '(A,I5)') 'Number of unique G-vectors: ', mpdata%num_gpts()
176 6 : WRITE(oUnit, *)
177 6 : WRITE(oUnit, '(3x,A)') 'IR Plane-wave basis with cutoff of gcutm (mpinp%g_cutoff/2*input%rkmax):'
178 54 : WRITE(oUnit, '(5x,A,I5)') 'Maximal number of G-vectors:', maxval(mpdata%n_g)
179 : END if
180 12 : end subroutine mpdata_gen_gvec
181 :
182 200 : subroutine mpdata_check_orthonormality(mpdata, atoms, mpi, l, itype, gridf)
183 :
184 : use m_judft
185 : use m_types_setup
186 : use m_types_mpi
187 : USE m_constants
188 : use m_intgrf, only: intgrf
189 :
190 : implicit none
191 :
192 : class(t_mpdata) :: mpdata
193 : type(t_atoms), intent(in) :: atoms
194 : type(t_mpi), intent(in) :: mpi
195 : integer, intent(in) :: itype, l
196 : real, intent(in) :: gridf(:, :)
197 :
198 200 : real, allocatable :: olap(:, :)
199 : integer :: i, n_radbasfn, err_loc(2)
200 :
201 : CHARACTER, PARAMETER :: lchar(0:38) = ['s', 'p', 'd', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', &
202 : 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', &
203 : 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x']
204 :
205 200 : call timestart("check mpdata orthonormality")
206 :
207 : ! calculate overlap matrix
208 200 : call mpdata%calc_olap_radbasfn(atoms, l, itype, gridf, olap)
209 :
210 : !subtract identity-matrix
211 1720 : do i = 1, size(olap, dim=1)
212 1720 : olap(i, i) = olap(i, i) - 1.0
213 : enddo
214 :
215 : ! check if (olap - identity) is zero-matrix
216 13608 : if(norm2(olap) > 1e-6) then
217 0 : if(mpi%irank == 0) THEN
218 0 : err_loc = maxloc(abs(olap))
219 : WRITE(*, *) 'mixedbasis: Bad orthonormality of ' &
220 0 : //lchar(l)//'-mixet product basis. Increase tolerance.'
221 0 : write(*, *) "itype =", itype, "l =", l
222 0 : WRITE(*, *) 'Deviation of', &
223 0 : maxval(abs(olap)), ' occurred for (', &
224 0 : err_loc(1), ',', err_loc(2), ')-overlap.'
225 : endif
226 : CALL judft_error("Bad orthonormality of mpdata", &
227 : hint='Increase tolerance', &
228 0 : calledby='mixedbasis%check_orthonormality')
229 : endif
230 :
231 200 : if(mpi%irank == 0) THEN
232 100 : n_radbasfn = mpdata%num_radbasfn(l, itype)
233 : WRITE(oUnit, '(6X,A,I4,'' ('',F22.18,'' )'')') &
234 6804 : lchar(l)//':', n_radbasfn, norm2(olap)/n_radbasfn
235 : END if
236 200 : call timestop("check mpdata orthonormality")
237 200 : end subroutine mpdata_check_orthonormality
238 :
239 12 : subroutine mpdata_check_radbasfn(mpdata, atoms, hybinp)
240 : use m_judft
241 : use m_types_hybinp
242 : use m_types_setup
243 : implicit none
244 : class(t_mpdata), intent(in) :: mpdata
245 : type(t_atoms), intent(in) :: atoms
246 : type(t_hybinp), intent(in) :: hybinp
247 :
248 : integer :: itype
249 :
250 32 : do itype = 1, atoms%ntype
251 132 : if(ANY(mpdata%num_radbasfn(0:hybinp%lcutm1(itype), itype) == 0)) THEN
252 0 : call judft_error('any mpdata%num_radbasfn eq 0', calledby='mixedbasis')
253 : endif
254 : enddo
255 12 : end subroutine mpdata_check_radbasfn
256 :
257 300 : SUBROUTINE mpdata_calc_olap_radbasfn(mpdata, atoms, l, itype, gridf, olap)
258 : USE ieee_arithmetic
259 : use m_intgrf, only: intgrf
260 : use m_types_setup
261 : use m_judft
262 : use m_types_fleurinput_base, only: REAL_NOT_INITALIZED
263 :
264 : implicit NONE
265 : class(t_mpdata), intent(in) :: mpdata
266 : type(t_atoms), intent(in) :: atoms
267 : integer, intent(in) :: l, itype
268 : real, intent(in) :: gridf(:, :)
269 : real, intent(inout), allocatable :: olap(:, :)
270 :
271 : integer :: n1, n2, n_radbasfn
272 :
273 300 : call timestart("calc mpdata overlap")
274 :
275 300 : n_radbasfn = mpdata%num_radbasfn(l, itype)
276 300 : if(allocated(olap)) then
277 88 : if(any(shape(olap) /= n_radbasfn)) then
278 88 : deallocate(olap)
279 : endif
280 : endif
281 300 : if(.not. allocated(olap)) allocate(olap(n_radbasfn, n_radbasfn), &
282 92960 : source=REAL_NOT_INITALIZED)
283 :
284 4316 : do n2 = 1, n_radbasfn
285 50196 : do n1 = 1, n2
286 : olap(n1, n2) = intgrf(mpdata%radbasfn_mt(:, n1, l, itype)*mpdata%radbasfn_mt(:, n2, l, itype), &
287 36131800 : atoms, itype, gridf)
288 : if(ieee_is_nan(olap(n1, n2))) then
289 : write(*, *) "nan at", n1, n2
290 : endif
291 49896 : olap(n2, n1) = olap(n1, n2)
292 : END do
293 : END do
294 : if(any(ieee_is_nan(olap))) call juDFT_error("Mixed-product basis olap is nan")
295 300 : call timestop("calc mpdata overlap")
296 300 : end subroutine mpdata_calc_olap_radbasfn
297 :
298 100 : subroutine mpdata_filter_radbasfn(mpdata, mpinp, l, itype, n_radbasfn, eig, eigv)
299 : ! Get rid of linear dependencies (eigenvalue <= mpdata%linear_dep_tol)
300 300 : use m_judft
301 : use m_types_mpinp
302 : implicit none
303 : class(t_mpdata), intent(inout) :: mpdata
304 : type(t_mpinp), intent(in) :: mpinp
305 : integer, intent(in) :: l, itype, n_radbasfn
306 : real, intent(inout) :: eig(:), eigv(:, :)
307 :
308 : integer :: num_radbasfn, i_bas
309 100 : integer, allocatable :: remaining_basfn(:)
310 :
311 100 : call timestart("filer mpdata")
312 2796 : allocate(remaining_basfn(n_radbasfn), source=1)
313 100 : num_radbasfn = 0
314 :
315 2596 : do i_bas = 1, mpdata%num_radbasfn(l, itype)
316 2596 : if(eig(i_bas) > mpinp%linear_dep_tol) THEN
317 740 : num_radbasfn = num_radbasfn + 1
318 740 : remaining_basfn(num_radbasfn) = i_bas
319 : END if
320 : END do
321 :
322 100 : mpdata%num_radbasfn(l, itype) = num_radbasfn
323 5292 : eig = eig(remaining_basfn)
324 156904 : eigv(:, :) = eigv(:, remaining_basfn)
325 100 : call timestop("filer mpdata")
326 100 : end subroutine mpdata_filter_radbasfn
327 :
328 100 : subroutine mpdata_diagonialize_olap(olap, eig_val, eig_vec)
329 : use m_judft
330 : use m_types_fleurinput_base, only: REAL_NOT_INITALIZED
331 : implicit NONE
332 : real, intent(in) :: olap(:, :)
333 : real, allocatable :: eig_val(:), eig_vec(:, :)
334 :
335 : integer :: n, size_iwork(1), info
336 : real :: size_work(1)
337 100 : integer, allocatable :: iwork(:)
338 100 : real, allocatable :: work(:)
339 :
340 100 : call timestart("diagonalize overlap")
341 100 : if(size(olap, dim=1) /= size(olap, dim=2)) then
342 0 : call juDFT_error("only square matrices can be diagonalized")
343 : endif
344 :
345 100 : n = size(olap, dim=1)
346 :
347 100 : if(allocated(eig_val)) then
348 88 : if(size(eig_val) /= n) deallocate(eig_val)
349 : endif
350 300 : if(.not. allocated(eig_val)) allocate(eig_val(n))
351 2596 : eig_val = REAL_NOT_INITALIZED
352 :
353 78552 : eig_vec = olap
354 : ! get sizes of work arrays
355 : call dsyevd('V', 'U', n, eig_vec, n, eig_val, &
356 100 : size_work, -1, size_iwork, -1, info)
357 100 : if(info /= 0) call juDFT_error("diagonalization for size failed")
358 :
359 300 : allocate(work(int(size_work(1))))
360 300 : allocate(iwork(size_iwork(1)))
361 :
362 : call dsyevd('V', 'U', n, eig_vec, n, eig_val, &
363 100 : work, int(size_work(1)), iwork, size_iwork(1), info)
364 100 : if(info /= 0) call juDFT_error("diagonalization failed")
365 100 : call timestop("diagonalize overlap")
366 100 : end subroutine mpdata_diagonialize_olap
367 :
368 100 : subroutine mpdata_trafo_to_orthonorm_bas(mpdata, full_n_radbasfn, n_grid_pt, l, itype, eig, eigv)
369 : use m_judft
370 : implicit NONE
371 : class(t_mpdata), intent(inout) :: mpdata
372 : integer, intent(in) :: full_n_radbasfn, n_grid_pt, l, itype
373 : real, intent(in) :: eig(:), eigv(:, :)
374 :
375 : integer :: nn, i
376 :
377 100 : call timestart("transform to reduced mpdata")
378 : ! reduced number of basis functions
379 100 : nn = mpdata%num_radbasfn(l, itype)
380 :
381 80000 : do i = 1, n_grid_pt
382 : mpdata%radbasfn_mt(i, 1:nn, l, itype) &
383 16906296 : = MATMUL(mpdata%radbasfn_mt(i, 1:full_n_radbasfn, l, itype), eigv(:, 1:nn))/SQRT(eig(:nn))
384 : END do
385 100 : call timestop("transform to reduced mpdata")
386 100 : end subroutine mpdata_trafo_to_orthonorm_bas
387 :
388 20 : subroutine mpdata_add_l0_fun(mpdata, atoms, hybinp, n_grid_pt, l, itype, gridf)
389 : use m_types_setup
390 : use m_types_hybinp
391 : use m_intgrf, only: intgrf
392 : use m_judft
393 : implicit none
394 : class(t_mpdata), intent(inout) :: mpdata
395 : type(t_atoms), intent(in) :: atoms
396 : type(t_hybinp), intent(in) :: hybinp
397 : integer, intent(in) :: n_grid_pt, l, itype
398 : real, intent(in) :: gridf(:, :)
399 :
400 : integer :: i, j, nn
401 20 : real, allocatable :: basmhlp(:, :, :, :)
402 : real :: norm
403 :
404 20 : call timestart("add l0 to mpdata")
405 20 : nn = mpdata%num_radbasfn(l, itype)
406 20 : if(l == 0) THEN
407 :
408 : ! Check if radbasfn_mt must be reallocated
409 20 : if(nn + 1 > SIZE(mpdata%radbasfn_mt, 2)) THEN
410 0 : allocate(basmhlp(atoms%jmtd, nn + 1, 0:maxval(hybinp%lcutm1), atoms%ntype))
411 0 : basmhlp(:, 1:nn, :, :) = mpdata%radbasfn_mt
412 0 : deallocate(mpdata%radbasfn_mt)
413 0 : allocate(mpdata%radbasfn_mt(atoms%jmtd, nn + 1, 0:maxval(hybinp%lcutm1), atoms%ntype))
414 0 : mpdata%radbasfn_mt(:, 1:nn, :, :) = basmhlp(:, 1:nn, :, :)
415 0 : deallocate(basmhlp)
416 : END if
417 :
418 : ! add l = 0 function
419 : mpdata%radbasfn_mt(:n_grid_pt, nn + 1, 0, itype) &
420 16000 : = atoms%rmsh(:n_grid_pt, itype)/SQRT(atoms%rmsh(n_grid_pt, itype)**3/3)
421 :
422 : ! perform gram-schmidt orthonormalization
423 180 : do i = nn, 1, -1
424 884 : do j = i + 1, nn + 1
425 : mpdata%radbasfn_mt(:n_grid_pt, i, 0, itype) &
426 : = mpdata%radbasfn_mt(:n_grid_pt, i, 0, itype) &
427 : - intgrf( &
428 : mpdata%radbasfn_mt(:n_grid_pt, i, 0, itype) &
429 : *mpdata%radbasfn_mt(:n_grid_pt, j, 0, itype), &
430 : atoms, itype, gridf) &
431 1152128 : *mpdata%radbasfn_mt(:n_grid_pt, j, 0, itype)
432 :
433 : END do
434 :
435 : ! renormalize
436 127784 : norm = SQRT(intgrf(mpdata%radbasfn_mt(:n_grid_pt, i, 0, itype)**2, atoms, itype, gridf))
437 : mpdata%radbasfn_mt(:n_grid_pt, i, 0, itype) &
438 127644 : = mpdata%radbasfn_mt(:n_grid_pt, i, 0, itype)/norm
439 : END do
440 20 : nn = nn + 1
441 20 : mpdata%num_radbasfn(l, itype) = nn
442 : END if
443 20 : call timestop("add l0 to mpdata")
444 20 : end subroutine mpdata_add_l0_fun
445 :
446 12 : subroutine mpdata_reduce_linear_dep(mpdata, mpinp, atoms, mpi, hybinp, gridf, iterHF)
447 : use m_types_setup
448 : use m_types_hybinp
449 : use m_types_mpi
450 : use m_judft
451 : use m_types_mpinp
452 : implicit none
453 : class(t_mpdata) :: mpdata
454 : type(t_mpinp), intent(in) :: mpinp
455 : type(t_atoms), intent(in) :: atoms
456 : type(t_mpi), intent(in) :: mpi
457 : type(t_hybinp), intent(in) :: hybinp
458 : integer, intent(in) :: iterHF
459 :
460 12 : real, allocatable :: olap(:, :), eig(:), eigv(:, :)
461 : real :: gridf(:, :)
462 : integer :: full_n_radbasfn, n_grid_pt, l, itype
463 :
464 : ! In order to get rid of the linear dependencies in the
465 : ! radial functions radbasfn_mt belonging to fixed l and itype
466 : ! the overlap matrix is diagonalized and those eigenvectors
467 : ! with a eigenvalue greater then mpdata%linear_dep_tol are retained
468 :
469 12 : call timestart("reduce lin. dep. mpdata")
470 :
471 32 : do itype = 1, atoms%ntype
472 20 : if (hybinp%lcutm1(itype) <= 0) call judft_error("lcutm1 <= 0 isn't allowed for hybrid calculations")
473 132 : DO l = 0, hybinp%lcutm1(itype)
474 100 : full_n_radbasfn = mpdata%num_radbasfn(l, itype)
475 100 : n_grid_pt = atoms%jri(itype)
476 :
477 : ! Calculate overlap
478 100 : call mpdata%calc_olap_radbasfn(atoms, l, itype, gridf, olap)
479 :
480 : ! Diagonalize
481 100 : call mpdata_diagonialize_olap(olap, eig, eigv)
482 :
483 100 : call mpdata%filter_radbasfn(mpinp, l, itype, full_n_radbasfn, eig, eigv)
484 :
485 100 : call mpdata%trafo_to_orthonorm_bas(full_n_radbasfn, n_grid_pt, l, itype, eig, eigv)
486 :
487 : ! Add constant function to l=0 basis and then do a Gram-Schmidt orthonormalization
488 100 : if(l == 0) then
489 20 : call mpdata%add_l0_fun(atoms, hybinp, n_grid_pt, l, itype, gridf)
490 : endif
491 :
492 : ! Check orthonormality of mpdatauct basis
493 120 : call mpdata%check_orthonormality(atoms, mpi, l, itype, gridf)
494 : enddo
495 : enddo
496 :
497 12 : deallocate(olap, eigv, eig)
498 12 : call timestop("reduce lin. dep. mpdata")
499 12 : end subroutine
500 :
501 12 : subroutine mpdata_normalize(mpdata, atoms, hybinp, gridf)
502 : use m_intgrf, only: intgrf
503 : use m_types_hybinp
504 : use m_types_setup
505 : use m_judft
506 : implicit NONE
507 :
508 : class(t_mpdata), intent(inout):: mpdata
509 : type(t_atoms), intent(in) :: atoms
510 : type(t_hybinp), intent(in) :: hybinp
511 : real, intent(in) :: gridf(:, :)
512 :
513 : integer :: l, i_basfn, itype
514 : real :: norm
515 :
516 32 : do itype = 1, atoms%ntype
517 132 : do l = 0, hybinp%lcutm1(itype)
518 2616 : do i_basfn = 1, mpdata%num_radbasfn(l, itype)
519 : norm = SQRT( &
520 : intgrf(mpdata%radbasfn_mt(:, i_basfn, l, itype)**2, &
521 : atoms, itype, gridf) &
522 1998176 : )
523 : mpdata%radbasfn_mt(:atoms%jri(itype), i_basfn, l, itype) &
524 1969572 : = mpdata%radbasfn_mt(:atoms%jri(itype), i_basfn, l, itype)/norm
525 : end do
526 : end do
527 : end do
528 12 : end subroutine mpdata_normalize
529 :
530 12 : subroutine mpdata_init(mpdata, hybinp, hybdat, atoms)
531 : use m_types_setup
532 : use m_types_hybinp
533 : use m_types_hybdat
534 : use m_judft
535 : implicit none
536 : class(t_mpdata) :: mpdata
537 : type(t_hybinp), intent(in) :: hybinp
538 : type(t_hybdat), intent(in) :: hybdat
539 : type(t_atoms), intent(in) :: atoms
540 :
541 : integer :: ok
542 :
543 :
544 12 : call mpdata%set_max_indx_p_1(atoms, hybinp)
545 :
546 12 : if(.not. allocated(mpdata%l1)) then
547 : allocate(mpdata%l1(mpdata%max_indx_p_1, 0:maxval(hybinp%lcutm1), atoms%ntype),&
548 40 : stat = ok)
549 6 : if(ok /= 0) call judft_error('mpdata_init: failure allocation mpdata%l1')
550 :
551 30 : allocate(mpdata%l2, mold=mpdata%l1, stat=ok)
552 6 : if(ok /= 0) call judft_error('mpdata_init: failure allocation mpdata%l2')
553 :
554 30 : allocate(mpdata%n1, mold=mpdata%l1, stat=ok)
555 6 : if(ok /= 0) call judft_error('mpdata_init: failure allocation mpdata%n1')
556 :
557 30 : allocate(mpdata%n2, mold=mpdata%l1, stat=ok)
558 6 : if(ok /= 0) call judft_error('mpdata_init: failure allocation mpdata%n2')
559 :
560 50406 : mpdata%l1 = -1; mpdata%l2 = -1; mpdata%n1 = -1; mpdata%n2 = -1
561 : endif
562 12 : end subroutine mpdata_init
563 :
564 12 : subroutine set_max_indx_p_1(mpdata, atoms, hybinp)
565 : use m_types_atoms
566 : use m_types_hybinp
567 : implicit none
568 : class(t_mpdata) :: mpdata
569 : type(t_atoms), intent(in) :: atoms
570 : type(t_hybinp), intent(in) :: hybinp
571 :
572 : integer :: itype, l, M, l1, l2, n1, n2
573 :
574 12 : mpdata%max_indx_p_1 = 0
575 :
576 32 : DO itype = 1, atoms%ntype
577 132 : DO l = 0, hybinp%lcutm1(itype)
578 100 : M = 0
579 1080 : DO l1 = 0, atoms%lmax(itype)
580 10780 : DO l2 = 0, atoms%lmax(itype)
581 10680 : IF(l >= ABS(l1 - l2) .AND. l <= l1 + l2) THEN
582 11396 : DO n1 = 1, mpdata%num_radfun_per_l(l1, itype)
583 27468 : DO n2 = 1, mpdata%num_radfun_per_l(l2, itype)
584 23768 : M = M + 1
585 : END DO
586 : END DO
587 : END IF
588 : END DO
589 : END DO
590 120 : mpdata%max_indx_p_1 = MAX(mpdata%max_indx_p_1, M)
591 : END DO
592 : END DO
593 12 : end subroutine set_max_indx_p_1
594 :
595 0 : subroutine mpdata_free(mpdata)
596 : implicit NONE
597 : class(t_mpdata) :: mpdata
598 :
599 0 : if(allocated(mpdata%l1)) deallocate(mpdata%l1)
600 0 : if(allocated(mpdata%l2)) deallocate(mpdata%l2)
601 0 : if(allocated(mpdata%n1)) deallocate(mpdata%n1)
602 0 : if(allocated(mpdata%n2)) deallocate(mpdata%n2)
603 0 : end subroutine mpdata_free
604 :
605 0 : subroutine mpdata_set_nl(mpdata, n, l, itype, n1, l1, n2, l2)
606 : implicit NONE
607 : class(t_mpdata) :: mpdata
608 : integer, intent(in) :: n, l, itype
609 : integer, intent(out) :: n1, l1, n2, l2
610 :
611 0 : l1 = mpdata%l1(n, l, itype) !
612 0 : l2 = mpdata%l2(n, l, itype) ! current basis-function mpdatauct
613 0 : n1 = mpdata%n1(n, l, itype) ! = bas(:,n1,l1,itype)*bas(:,n2,l2,itype) = b1*b2
614 0 : n2 = mpdata%n2(n, l, itype) !
615 :
616 0 : end subroutine mpdata_set_nl
617 :
618 12 : subroutine set_num_radfun_per_l_mpdata(mpdata, atoms)
619 : use m_types_setup
620 : implicit NONE
621 : class(t_mpdata) :: mpdata
622 : type(t_atoms) :: atoms
623 : integer :: itype, ilo
624 :
625 12 : if(.not. allocated(mpdata%num_radfun_per_l)) then
626 0 : allocate(mpdata%num_radfun_per_l(0:atoms%lmaxd, atoms%ntype))
627 : endif
628 :
629 : ! always two there are: u and u_dot
630 228 : mpdata%num_radfun_per_l = 2
631 32 : DO itype = 1, atoms%ntype
632 56 : DO ilo = 1, atoms%nlo(itype)
633 : mpdata%num_radfun_per_l(atoms%llo(ilo, itype), itype) &
634 44 : = mpdata%num_radfun_per_l(atoms%llo(ilo, itype), itype) + 1
635 : END DO
636 : END DO
637 12 : end subroutine set_num_radfun_per_l_mpdata
638 79900 : end module m_types_mpdata
|
