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 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
8 : ! This subroutine generates the mixed basis set used to evaluate the !
9 : ! exchange term in HF/hybinp functional calculations or EXX !
10 : ! calculations. In the latter case a second mixed basis set is setup !
11 : ! for the OEP integral equation. !
12 : ! In all cases the mixed basis consists of IR plane waves !
13 : ! !
14 : ! IR: !
15 : ! M_{\vec{k},\vec{G}} = 1/\sqrt{V} \exp{i(\vec{k}+\vec{G})} !
16 : ! !
17 : ! which are zero in the MT spheres and radial functions times !
18 : ! spherical harmonics in the MT spheres !
19 : ! !
20 : ! MT: !
21 : ! a a !
22 : ! M = U Y !
23 : ! PL PL LM !
24 : ! !
25 : ! where a a a !
26 : ! U = u u !
27 : ! PL pl p'l' !
28 : ! !
29 : ! and L \in {|l-l'|,...,l+l'} !
30 : ! !
31 : ! and P counts the different combinations of !
32 : ! pl,p'l' which contribute to L !
33 : ! !
34 : ! M.Betzinger (09/07) !
35 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
36 :
37 : MODULE m_mixedbasis
38 :
39 : CONTAINS
40 :
41 12 : SUBROUTINE mixedbasis(atoms, kpts, input, cell, xcpot, mpinp, mpdata, hybinp, hybdat,&
42 : enpara, fmpi, v, iterHF)
43 :
44 : USE m_judft
45 : USE m_types
46 : USE m_constants
47 : USE m_loddop, ONLY: loddop
48 : USE m_intgrf, ONLY: intgrf_init, intgrf
49 : USE m_hybrid_core
50 : USE m_wrapper
51 : USE m_eig66_io
52 :
53 : IMPLICIT NONE
54 :
55 : TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
56 : TYPE(t_mpi), INTENT(IN) :: fmpi
57 : TYPE(t_mpdata), intent(inout) :: mpdata
58 : TYPE(t_mpinp), intent(in) :: mpinp
59 : TYPE(t_hybinp), INTENT(IN) :: hybinp
60 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat
61 : TYPE(t_enpara), INTENT(IN) :: enpara
62 : TYPE(t_input), INTENT(IN) :: input
63 : TYPE(t_cell), INTENT(IN) :: cell
64 : TYPE(t_kpts), INTENT(IN) :: kpts
65 : TYPE(t_atoms), INTENT(IN) :: atoms
66 : TYPE(t_potden), INTENT(IN) :: v
67 :
68 : integer, intent(in) :: iterHF
69 :
70 : ! local type variables
71 12 : TYPE(t_usdus) :: usdus
72 :
73 : ! local scalars
74 : INTEGER :: jspin, itype, l1, l2, l, n_radbasfn, full_n_radbasfn, n1, n2
75 : INTEGER :: i_basfn, i, n_grid_pt,j
76 : REAL :: rdum, rdum1, max_momentum, momentum
77 :
78 : ! - local arrays -
79 :
80 12 : REAL :: bashlp(atoms%jmtd)
81 :
82 12 : REAL, ALLOCATABLE :: bas1(:,:,:,:,:), bas2(:,:,:,:,:)
83 12 : REAL, ALLOCATABLE :: basmhlp(:,:,:,:)
84 12 : REAL, ALLOCATABLE :: gridf(:,:), vr0(:,:,:)
85 :
86 12 : LOGICAL, ALLOCATABLE :: selecmat(:,:,:,:)
87 12 : LOGICAL, ALLOCATABLE :: seleco(:,:), selecu(:,:)
88 :
89 : CHARACTER, PARAMETER :: lchar(0:38) = (/'s', 'p', 'd', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', &
90 : 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', &
91 : 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x'/)
92 :
93 :
94 12 : IF (fmpi%irank == 0) WRITE (oUnit, '(//A,I2,A)') '### subroutine: mixedbasis ###'
95 :
96 12 : IF (xcpot%is_name("exx")) CALL judft_error("EXX is not implemented in this version", calledby='mixedbasis')
97 :
98 : ! Deallocate arrays which might have been allocated in a previous run of this subroutine
99 12 : IF (ALLOCATED(mpdata%n_g)) deallocate(mpdata%n_g)
100 12 : IF (ALLOCATED(mpdata%num_radbasfn)) deallocate(mpdata%num_radbasfn)
101 12 : IF (ALLOCATED(mpdata%gptm_ptr)) deallocate(mpdata%gptm_ptr)
102 12 : IF (ALLOCATED(mpdata%g)) deallocate(mpdata%g)
103 12 : IF (ALLOCATED(mpdata%radbasfn_mt)) deallocate(mpdata%radbasfn_mt)
104 :
105 12 : CALL usdus%init(atoms, input%jspins)
106 :
107 12 : if(.not. allocated(mpdata%num_radfun_per_l)) &
108 132 : allocate(mpdata%num_radfun_per_l(0:atoms%lmaxd, atoms%ntype), source=-1)
109 :
110 12 : call mpdata%set_num_radfun_per_l(atoms)
111 12 : call mpdata%init(hybinp, hybdat, atoms)
112 :
113 : ! initialize gridf for radial integration
114 12 : CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
115 :
116 19404 : allocate(vr0(atoms%jmtd, atoms%ntype, input%jspins), source=0.0)
117 :
118 19356 : vr0(:,:,:) = v%mt(:,0, :,:)
119 :
120 : ! calculate radial basisfunctions u and u' with
121 : ! the spherical part of the potential vr0 and store them in
122 : ! bas1 = large component ,bas2 = small component
123 :
124 12 : call gen_bas_fun(atoms, enpara, gridf, input, mpdata, fmpi, vr0, usdus, bas1, bas2)
125 :
126 : ! - - - - - - SETUP OF THE MIXED BASIS IN THE IR - - - - - - -
127 :
128 : ! construct G-vectors with cutoff smaller than gcutm
129 12 : call mpdata%gen_gvec(mpinp, cell, kpts, fmpi)
130 :
131 : ! - - - - - - - - Set up MT product basis for the non-local exchange potential - - - - - - - - - -
132 :
133 12 : IF (fmpi%irank == 0) THEN
134 6 : WRITE (oUnit, '(A)') 'MT product basis for non-local exchange potential:'
135 6 : WRITE (oUnit, '(A)') 'Reduction due to overlap (quality of orthonormality, should be < 1.0E-06)'
136 : END IF
137 :
138 188 : allocate(mpdata%num_radbasfn(0:maxval(hybinp%lcutm1), atoms%ntype), source=0)
139 264 : allocate(seleco(maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd))
140 264 : allocate(selecu(maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd))
141 :
142 : ! determine maximal indices of (radial) mixed-basis functions (->num_radbasfn)
143 : ! (will be reduced later-on due to overlap)
144 32 : DO itype = 1, atoms%ntype
145 804 : seleco = .FALSE.
146 804 : selecu = .FALSE.
147 120 : seleco(1, 0:hybinp%select1(1, itype)) = .TRUE.
148 120 : selecu(1, 0:hybinp%select1(3, itype)) = .TRUE.
149 40 : seleco(2, 0:hybinp%select1(2, itype)) = .TRUE.
150 80 : selecu(2, 0:hybinp%select1(4, itype)) = .TRUE.
151 :
152 : ! include local orbitals
153 412 : IF (maxval(mpdata%num_radfun_per_l) >= 3) THEN
154 412 : seleco(3:,:) = .TRUE.
155 412 : selecu(3:,:) = .TRUE.
156 : END IF
157 :
158 132 : DO l = 0, hybinp%lcutm1(itype)
159 100 : n_radbasfn = 0
160 :
161 : !
162 : ! valence * valence
163 : !
164 100 : if(.not. allocated(selecmat)) then
165 0 : allocate(selecmat(maxval(mpdata%num_radfun_per_l), &
166 : 0:atoms%lmaxd, &
167 : maxval(mpdata%num_radfun_per_l), &
168 288 : 0:atoms%lmaxd))
169 : endif
170 4020 : selecmat = calc_selecmat(atoms, mpdata, seleco, selecu)
171 :
172 1080 : DO l1 = 0, atoms%lmax(itype)
173 10780 : DO l2 = 0, atoms%lmax(itype)
174 10680 : IF (l >= ABS(l1 - l2) .AND. l <= l1 + l2) THEN
175 11396 : DO n1 = 1, mpdata%num_radfun_per_l(l1, itype)
176 27468 : DO n2 = 1, mpdata%num_radfun_per_l(l2, itype)
177 23768 : IF (selecmat(n1, l1, n2, l2)) THEN
178 2056 : n_radbasfn = n_radbasfn + 1
179 2056 : selecmat(n2, l2, n1, l1) = .FALSE. ! prevent double counting of products (a*b = b*a)
180 : END IF
181 : END DO
182 : END DO
183 : END IF
184 : END DO
185 : END DO
186 100 : IF (n_radbasfn == 0 .AND. fmpi%irank == 0) &
187 : WRITE (oUnit, '(A)') 'mixedbasis: Warning! No basis-function product of '//lchar(l)// &
188 0 : '-angular momentum defined.'
189 120 : mpdata%num_radbasfn(l, itype) = n_radbasfn*input%jspins
190 : END DO
191 : END DO
192 :
193 0 : allocate(mpdata%radbasfn_mt(atoms%jmtd,&
194 : maxval(mpdata%num_radbasfn), &
195 : 0:maxval(hybinp%lcutm1), &
196 2735932 : atoms%ntype), source=0.0)
197 :
198 : ! Define product bases and reduce them according to overlap
199 :
200 32 : DO itype = 1, atoms%ntype
201 804 : seleco = .FALSE.
202 804 : selecu = .FALSE.
203 120 : seleco(1, 0:hybinp%select1(1, itype)) = .TRUE.
204 120 : selecu(1, 0:hybinp%select1(3, itype)) = .TRUE.
205 40 : seleco(2, 0:hybinp%select1(2, itype)) = .TRUE.
206 80 : selecu(2, 0:hybinp%select1(4, itype)) = .TRUE.
207 : ! include lo's
208 412 : IF (maxval(mpdata%num_radfun_per_l) >= 3) THEN
209 412 : seleco(3:,:) = .TRUE.
210 412 : selecu(3:,:) = .TRUE.
211 : END IF
212 :
213 20 : n_grid_pt = atoms%jri(itype)
214 120 : DO l = 0, hybinp%lcutm1(itype)
215 100 : full_n_radbasfn = mpdata%num_radbasfn(l, itype)
216 : ! allow for zero product-basis functions for
217 : ! current l-quantum number
218 100 : IF (n_radbasfn == 0) THEN
219 0 : IF (fmpi%irank == 0) WRITE (oUnit, '(6X,A,'': 0 -> 0'')') lchar(l)
220 : CYCLE
221 : END IF
222 :
223 : ! set up the overlap matrix
224 100 : i_basfn = 0
225 :
226 : ! valence*valence
227 4020 : selecmat = calc_selecmat(atoms, mpdata,seleco, selecu)
228 :
229 1080 : DO l1 = 0, atoms%lmax(itype)
230 10780 : DO l2 = 0, atoms%lmax(itype)
231 10680 : IF (l >= ABS(l1 - l2) .AND. l <= l1 + l2) THEN
232 11396 : DO n1 = 1, mpdata%num_radfun_per_l(l1, itype)
233 27468 : DO n2 = 1, mpdata%num_radfun_per_l(l2, itype)
234 :
235 23768 : IF (selecmat(n1, l1, n2, l2)) THEN
236 4552 : DO jspin = 1, input%jspins
237 2496 : i_basfn = i_basfn + 1
238 2496 : IF (i_basfn > full_n_radbasfn) call judft_error('got too many product functions', hint='This is a BUG, please report', calledby='mixedbasis')
239 :
240 : mpdata%radbasfn_mt(:n_grid_pt, i_basfn, l, itype) &
241 : = ( bas1(:n_grid_pt, n1, l1, itype, jspin) &
242 : * bas1(:n_grid_pt, n2, l2, itype, jspin) &
243 : + bas2(:n_grid_pt, n1, l1, itype, jspin) &
244 : * bas2(:n_grid_pt, n2, l2, itype, jspin) &
245 1971528 : ) / atoms%rmsh(:n_grid_pt, itype)
246 :
247 :
248 : END DO !jspin
249 : ! prevent double counting of products (a*b = b*a)
250 2056 : selecmat(n2, l2, n1, l1) = .FALSE.
251 : END IF
252 : END DO !n2
253 : END DO !n1
254 : ENDIF
255 : END DO !l2
256 : END DO !l1
257 :
258 :
259 120 : IF (i_basfn /= full_n_radbasfn) call judft_error('counting error for product functions', hint='This is a BUG, please report', calledby='mixedbasis')
260 :
261 :
262 : END DO !l
263 82 : IF (fmpi%irank == 0) WRITE (oUnit, '(6X,A,I7)') 'Total:', SUM(mpdata%num_radbasfn(0:hybinp%lcutm1(itype), itype))
264 : END DO ! itype
265 :
266 : !normalize radbasfn_mt
267 12 : call mpdata%normalize(atoms, hybinp, gridf)
268 12 : call mpdata%reduce_linear_dep(mpinp,atoms, fmpi, hybinp, gridf, iterHF)
269 :
270 212 : allocate(basmhlp(atoms%jmtd, maxval(mpdata%num_radbasfn), 0:maxval(hybinp%lcutm1), atoms%ntype))
271 : basmhlp(1:atoms%jmtd, 1:maxval(mpdata%num_radbasfn), 0:maxval(hybinp%lcutm1), 1:atoms%ntype) &
272 750352 : = mpdata%radbasfn_mt(1:atoms%jmtd, 1:maxval(mpdata%num_radbasfn), 0:maxval(hybinp%lcutm1), 1:atoms%ntype)
273 12 : deallocate(mpdata%radbasfn_mt)
274 212 : allocate(mpdata%radbasfn_mt(atoms%jmtd, maxval(mpdata%num_radbasfn), 0:maxval(hybinp%lcutm1), atoms%ntype))
275 750224 : mpdata%radbasfn_mt = basmhlp
276 :
277 12 : deallocate(basmhlp, seleco, selecu, selecmat)
278 :
279 : !
280 : ! now we build linear combinations of the radial functions
281 : ! such that they possess no moment except one radial function in each l-channel
282 : !
283 12 : IF (fmpi%irank == 0) THEN
284 : WRITE (oUnit, '(/,A,/,A)') 'Build linear combinations of radial '// &
285 6 : 'functions in each l-channel,', &
286 : 'such that they possess no multipolmoment'// &
287 12 : ' except the last function:'
288 :
289 6 : WRITE (oUnit, '(/,17x,A)') 'moment (quality of orthonormality)'
290 : END IF
291 :
292 32 : DO itype = 1, atoms%ntype
293 20 : n_grid_pt = atoms%jri(itype)
294 :
295 20 : IF (atoms%ntype > 1 .AND. fmpi%irank == 0) WRITE (oUnit, '(6X,A,I3)') 'Atom type', itype
296 :
297 120 : DO l = 0, hybinp%lcutm1(itype)
298 : ! determine radial function with the largest moment
299 : ! this function is used to build the linear combinations
300 100 : max_momentum = 0
301 860 : DO i = 1, mpdata%num_radbasfn(l, itype)
302 : momentum = intgrf(atoms%rmsh(:n_grid_pt, itype)**(l + 1)*mpdata%radbasfn_mt(:n_grid_pt, i, l, itype), &
303 606832 : atoms, itype, gridf)
304 860 : IF (ABS(momentum) > max_momentum) THEN
305 664 : n_radbasfn = i
306 664 : max_momentum = momentum
307 : END IF
308 : END DO
309 :
310 100 : j = 0
311 80000 : bashlp(:atoms%jri(itype)) = mpdata%radbasfn_mt(:atoms%jri(itype), n_radbasfn, l, itype)
312 860 : DO i = 1, mpdata%num_radbasfn(l, itype)
313 760 : IF (i == n_radbasfn) CYCLE
314 660 : j = j + 1
315 526172 : mpdata%radbasfn_mt(:atoms%jri(itype), j, l, itype) = mpdata%radbasfn_mt(:atoms%jri(itype), i, l, itype)
316 : END DO
317 80020 : mpdata%radbasfn_mt(:atoms%jri(itype), mpdata%num_radbasfn(l, itype), l, itype) = bashlp(:atoms%jri(itype))
318 : END DO
319 :
320 :
321 132 : DO l = 0, hybinp%lcutm1(itype)
322 100 : IF (fmpi%irank == 0) WRITE (oUnit, '(6X,A)') lchar(l)//':'
323 :
324 100 : IF (mpdata%num_radbasfn(l, itype) == 0) THEN
325 0 : IF (fmpi%irank == 0) WRITE (oUnit, '(6X,A,'': 0 -> '')') lchar(l)
326 : CYCLE
327 : END IF
328 :
329 100 : n_radbasfn = mpdata%num_radbasfn(l, itype)
330 760 : DO i = 1, n_radbasfn-1
331 : ! calculate moment of radial function i
332 : rdum1 = intgrf(atoms%rmsh(:n_grid_pt, itype)**(l + 1)*mpdata%radbasfn_mt(:n_grid_pt, i, l, itype), &
333 526732 : atoms, itype, gridf)
334 :
335 : rdum = intgrf(atoms%rmsh(:n_grid_pt, itype)**(l + 1)*mpdata%radbasfn_mt(:n_grid_pt, n_radbasfn, l, itype), &
336 526072 : atoms, itype, gridf)
337 :
338 526072 : bashlp(:n_grid_pt) = mpdata%radbasfn_mt(:n_grid_pt, n_radbasfn, l, itype)
339 :
340 660 : IF (SQRT(rdum**2 + rdum1**2) <= 1E-06 .AND. fmpi%irank == 0) &
341 0 : WRITE (oUnit, *) 'Warning: Norm is smaller than 1E-06!'
342 :
343 : ! change function n_radbasfn such that n_radbasfn is orthogonal to i
344 : ! since the functions radbasfn_mt have been orthogonal on input
345 : ! the linear combination does not destroy the orthogonality to the residual functions
346 : mpdata%radbasfn_mt(:n_grid_pt, n_radbasfn, l, itype) = rdum/SQRT(rdum**2 + rdum1**2)*bashlp(:n_grid_pt) &
347 526072 : + rdum1/SQRT(rdum**2 + rdum1**2)*mpdata%radbasfn_mt(:n_grid_pt, i, l, itype)
348 :
349 : ! combine basis function i and n_radbasfn so that they possess no momemt
350 : mpdata%radbasfn_mt(:n_grid_pt, i, l, itype) = rdum1/SQRT(rdum**2 + rdum1**2)*bashlp(:n_grid_pt) &
351 526072 : - rdum/SQRT(rdum**2 + rdum1**2)*mpdata%radbasfn_mt(:n_grid_pt, i, l, itype)
352 :
353 : rdum1 = intgrf(atoms%rmsh(:n_grid_pt, itype)**(l + 1)*mpdata%radbasfn_mt(:n_grid_pt, i, l, itype), &
354 526072 : atoms, itype, gridf)
355 :
356 660 : IF (rdum1 > 1E-10) call judft_error('moment of radial function does not vanish', calledby='mixedbasis')
357 :
358 760 : IF (fmpi%irank == 0) WRITE (oUnit, '(6x,I4,'' -> '',F22.18)') i, rdum1
359 : END DO
360 120 : call mpdata%check_orthonormality(atoms, fmpi, l, itype, gridf)
361 : ENDDO
362 :
363 :
364 : END DO
365 :
366 12 : call mpdata%check_radbasfn(atoms, hybinp)
367 :
368 : !count basis functions
369 12 : hybdat%nbasp = 0
370 32 : DO itype = 1, atoms%ntype
371 52 : DO i = 1, atoms%neq(itype)
372 140 : DO l = 0, hybinp%lcutm1(itype)
373 120 : hybdat%nbasp = hybdat%nbasp + (2*l+1) * mpdata%num_radbasfn(l, itype)
374 : END DO
375 : END DO
376 : END DO
377 120 : hybdat%nbasm = hybdat%nbasp + mpdata%n_g
378 :
379 12 : hybdat%maxlmindx = 0
380 32 : do itype = 1,atoms%ntype
381 : hybdat%maxlmindx = max(hybdat%maxlmindx,&
382 : SUM([(mpdata%num_radfun_per_l(l, itype)*(2*l + 1), l=0, atoms%lmax(itype))])&
383 424 : )
384 : enddo
385 24 : END SUBROUTINE mixedbasis
386 :
387 12 : subroutine gen_bas_fun(atoms, enpara, gridf, input, mpdata, fmpi, vr0, usdus, bas1, bas2)
388 : use m_judft
389 : use m_types
390 : USE m_radfun, ONLY: radfun
391 : USE m_radflo, ONLY: radflo
392 : USE m_intgrf, ONLY: intgrf
393 : implicit NONE
394 : type(t_atoms), intent(in) :: atoms
395 : type(t_enpara), intent(in) :: enpara
396 : type(t_input), intent(in) :: input
397 : TYPE(t_mpdata), intent(in) :: mpdata
398 : type(t_mpi), intent(in) :: fmpi
399 : type(t_usdus), intent(inout) :: usdus
400 :
401 : REAL, ALLOCATABLE, INTENT(INOUT) :: bas1(:,:,:,:,:), bas2(:,:,:,:,:)
402 : REAL, intent(in) :: vr0(:,:,:), gridf(:,:)
403 :
404 12 : REAL :: u(atoms%jmtd, 2, 0:atoms%lmaxd)
405 12 : REAL :: du(atoms%jmtd, 2, 0:atoms%lmaxd)
406 12 : REAL :: flo(atoms%jmtd, 2, atoms%nlod)
407 :
408 12 : REAL :: uuilon(atoms%nlod, atoms%ntype)
409 12 : REAL :: duilon(atoms%nlod, atoms%ntype)
410 12 : REAL :: ulouilopn(atoms%nlod, atoms%nlod, atoms%ntype)
411 : REAL :: wronk, norm
412 :
413 : INTEGER :: itype, jspin, i, l, ilo, ok
414 : INTEGER :: n_grid_pt, noded, nodem
415 12 : INTEGER :: l_idx(0:atoms%lmaxd)
416 188896 : u = 0.0
417 188896 : du = 0.0
418 :
419 : ! this is 5-D array. it could cause Problems in bigger systems
420 : allocate(bas1(atoms%jmtd, &
421 : maxval(mpdata%num_radfun_per_l), &
422 : 0:atoms%lmaxd, &
423 : atoms%ntype, &
424 566876 : input%jspins), source=0.0, stat=ok)
425 12 : if(ok /= 0) call judft_error("Can't allocate bas1 array. Stat= " // int2str(ok))
426 :
427 566660 : allocate(bas2, source=bas1, stat=ok)
428 12 : if(ok /= 0) call judft_error("Can't allocate bas1 array. Stat= " // int2str(ok))
429 :
430 32 : DO itype = 1, atoms%ntype
431 20 : n_grid_pt = atoms%jri(itype) ! number of radial gridpoints
432 56 : DO jspin = 1, input%jspins
433 256 : DO l = 0, atoms%lmax(itype)
434 : CALL radfun(l, itype, jspin, enpara%el0(l, itype, jspin), vr0(:,itype, jspin), atoms, &
435 256 : u(:,:,l), du(:,:,l), usdus, nodem, noded, wronk)
436 : END DO
437 185096 : bas1(1:n_grid_pt, 1, 0:atoms%lmaxd, itype, jspin) = u(1:n_grid_pt, 1, 0:atoms%lmaxd)
438 185096 : bas2(1:n_grid_pt, 1, 0:atoms%lmaxd, itype, jspin) = u(1:n_grid_pt, 2, 0:atoms%lmaxd)
439 185096 : bas1(1:n_grid_pt, 2, 0:atoms%lmaxd, itype, jspin) = du(1:n_grid_pt, 1, 0:atoms%lmaxd)
440 185096 : bas2(1:n_grid_pt, 2, 0:atoms%lmaxd, itype, jspin) = du(1:n_grid_pt, 2, 0:atoms%lmaxd)
441 :
442 : ! generate radial functions for local orbitals
443 44 : IF (atoms%nlo(itype) >= 1) THEN
444 : CALL radflo(atoms, itype, jspin, enpara%ello0(1, 1, jspin), vr0(:,itype, jspin), &
445 20 : u, du, fmpi, usdus, uuilon, duilon, ulouilopn, flo)
446 :
447 208 : l_idx = 2
448 52 : DO ilo = 1, atoms%nlo(itype)
449 32 : l = atoms%llo(ilo, itype)
450 32 : l_idx(l) = l_idx(l) + 1
451 25312 : bas1(1:n_grid_pt, l_idx(l), atoms%llo(ilo, itype), itype, jspin) = flo(1:n_grid_pt, 1, ilo)
452 25332 : bas2(1:n_grid_pt, l_idx(l), atoms%llo(ilo, itype), itype, jspin) = flo(1:n_grid_pt, 2, ilo)
453 : END DO
454 : END IF
455 : END DO
456 : END DO
457 :
458 : ! the radial functions are normalized
459 28 : DO jspin = 1, input%jspins
460 52 : DO itype = 1, atoms%ntype
461 272 : DO l = 0, atoms%lmax(itype)
462 752 : DO i = 1, mpdata%num_radfun_per_l(l, itype)
463 : norm = sqrt(intgrf(bas1(:,i, l, itype, jspin)**2 + bas2(:,i, l, itype, jspin)**2, &
464 403344 : atoms, itype, gridf))
465 395456 : bas1(:atoms%jri(itype), i, l, itype, jspin) = bas1(:atoms%jri(itype), i, l, itype, jspin)/norm
466 395688 : bas2(:atoms%jri(itype), i, l, itype, jspin) = bas2(:atoms%jri(itype), i, l, itype, jspin)/norm
467 : END DO
468 : END DO
469 : END DO
470 : END DO
471 12 : end subroutine gen_bas_fun
472 :
473 8040 : function calc_selecmat(atoms,mpdata,seleco, selecu) result(selecmat)
474 : ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
475 : use m_types
476 : use m_judft
477 : implicit NONE
478 :
479 : type(t_atoms), intent(in) :: atoms
480 : TYPE(t_mpdata), intent(in) :: mpdata
481 : LOGICAL, intent(in) :: seleco(maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd)
482 : LOGICAL, intent(in) :: selecu(maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd)
483 : LOGICAL :: selecmat(maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, &
484 : maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd)
485 : integer :: n1, l1, n2, l2
486 :
487 : ! column-major means left-most index varies the fastest
488 2160 : do l2=0,atoms%lmaxd
489 49120 : do n2=1,maxval(mpdata%num_radfun_per_l)
490 66040 : do l1=0,atoms%lmaxd
491 1479840 : do n1=1,maxval(mpdata%num_radfun_per_l)
492 350960 : selecmat(n1,l1,n2,l2) = seleco(n1, l1) .AND. selecu(n2, l2)
493 : enddo
494 : enddo
495 : enddo
496 : enddo
497 200 : end function calc_selecmat
498 : END MODULE m_mixedbasis
|