mixedbasis.F90 27.8 KB
Newer Older
1
2
3
4
5
6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This subroutine generates the mixed basis set used to evaluate the  !
! exchange term in HF/hybrid functional calculations or EXX           !
! calculations. In the latter case a second mixed basis set is setup  !
! for the OEP integral equation.                                      !
! In all cases the mixed basis consists of IR plane waves             !
!                                                                     !
! IR:                                                                 !
!    M_{\vec{k},\vec{G}} = 1/\sqrt{V} \exp{i(\vec{k}+\vec{G})}        !
!                                                                     !
! which are zero in the MT spheres and radial functions times         !
! spherical harmonics in the MT spheres                               !
!                                                                     !
! MT:                                                                 !
!     a                a                                              !
!    M              = U   Y                                           !
!     PL               PL  LM                                         !
!                                                                     !
!            where     a    a  a                                      !
!                     U  = u  u                                       !
!                      PL   pl  p'l'                                  !
!                                                                     !
!               and    L \in {|l-l'|,...,l+l'}                        !
!                                                                     !
!               and    P counts the different combinations of         !
!                      pl,p'l' which contribute to L                  !
!                                                                     !
!                                               M.Betzinger (09/07)   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
36

Daniel Wortmann's avatar
Daniel Wortmann committed
37
38
MODULE m_mixedbasis

39
CONTAINS
40

Matthias Redies's avatar
Matthias Redies committed
41
   SUBROUTINE mixedbasis(atoms, kpts, input, cell, xcpot, mpbasis, hybrid, enpara, mpi, v)
Matthias Redies's avatar
Matthias Redies committed
42
43
44

      USE m_judft
      USE m_loddop, ONLY: loddop
Matthias Redies's avatar
Matthias Redies committed
45
      USE m_intgrf, ONLY: intgrf_init, intgrf
46
      use m_rorder, only: rorderpf
47
      USE m_hybrid_core
Matthias Redies's avatar
Matthias Redies committed
48
49
50
51
52
53
54
55
      USE m_wrapper
      USE m_eig66_io
      USE m_types

      IMPLICIT NONE

      TYPE(t_xcpot_inbuild), INTENT(IN)    :: xcpot
      TYPE(t_mpi), INTENT(IN)    :: mpi
Matthias Redies's avatar
Matthias Redies committed
56
      TYPE(t_mpbasis), intent(inout)  :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
57
58
59
60
61
62
63
64
65
66
67
68
69
      TYPE(t_hybrid), INTENT(INOUT) :: hybrid
      TYPE(t_enpara), INTENT(IN)    :: enpara
      TYPE(t_input), INTENT(IN)    :: input
      TYPE(t_cell), INTENT(IN)    :: cell
      TYPE(t_kpts), INTENT(IN)    :: kpts
      TYPE(t_atoms), INTENT(IN)    :: atoms
      TYPE(t_potden), INTENT(IN)    :: v


      ! local type variables
      TYPE(t_usdus)                   ::  usdus

      ! local scalars
70
71
      INTEGER                         ::  jspin, itype, l1, l2, l, n_radbasfn, n1, n2, nn
      INTEGER                         ::  m, nk, i_basfn, i, j, n_grid_pt
72
      REAL                            ::  rdum, rdum1, norm, max_momentum, momentum
Matthias Redies's avatar
Matthias Redies committed
73
74
75
76
77

      ! - local arrays -

      REAL                            ::  bashlp(atoms%jmtd)

Matthias Redies's avatar
Matthias Redies committed
78

Matthias Redies's avatar
Matthias Redies committed
79
80
81
82
      REAL, ALLOCATABLE               ::  olap(:,:), work(:), eig(:), eigv(:,:)
      REAL, ALLOCATABLE               ::  bas1(:,:,:,:,:), bas2(:,:,:,:,:)
      REAL, ALLOCATABLE               ::  basmhlp(:,:,:,:)
      REAL, ALLOCATABLE               ::  gridf(:,:), vr0(:,:,:)
Matthias Redies's avatar
Matthias Redies committed
83

Matthias Redies's avatar
Matthias Redies committed
84
85
      LOGICAL, ALLOCATABLE            ::  selecmat(:,:,:,:)
      LOGICAL, ALLOCATABLE            ::  seleco(:,:), selecu(:,:)
Matthias Redies's avatar
Matthias Redies committed
86
87
88
89
90
91
92
93
94
95
96

      CHARACTER, PARAMETER            :: lchar(0:38) = (/'s', 'p', 'd', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', &
                                                         'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', &
                                                         'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x'/)


      IF (mpi%irank == 0) WRITE (6, '(//A,I2,A)') '### subroutine: mixedbasis ###'

      IF (xcpot%is_name("exx")) CALL judft_error("EXX is not implemented in this version", calledby='mixedbasis')

      ! Deallocate arrays which might have been allocated in a previous run of this subroutine
Matthias Redies's avatar
Matthias Redies committed
97
      IF (ALLOCATED(mpbasis%ngptm)) deallocate(mpbasis%ngptm)
98
      IF (ALLOCATED(mpbasis%num_radbasfn)) deallocate(mpbasis%num_radbasfn)
Matthias Redies's avatar
Matthias Redies committed
99
      IF (ALLOCATED(mpbasis%gptm_ptr)) deallocate(mpbasis%gptm_ptr)
Matthias Redies's avatar
Matthias Redies committed
100
      IF (ALLOCATED(mpbasis%gptm)) deallocate(mpbasis%gptm)
Matthias Redies's avatar
Matthias Redies committed
101
      IF (ALLOCATED(mpbasis%radbasfn_mt)) deallocate(mpbasis%radbasfn_mt)
Matthias Redies's avatar
Matthias Redies committed
102
103

      CALL usdus%init(atoms, input%jspins)
104
      call hybrid%set_num_radfun_per_l(atoms)
Matthias Redies's avatar
Matthias Redies committed
105
106
107

      ! initialize gridf for radial integration
      CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
Daniel Wortmann's avatar
Daniel Wortmann committed
108

109
      allocate(vr0(atoms%jmtd, atoms%ntype, input%jspins), source=0.0)
Daniel Wortmann's avatar
Daniel Wortmann committed
110

Matthias Redies's avatar
Matthias Redies committed
111
      vr0(:,:,:) = v%mt(:,0, :,:)
112

Matthias Redies's avatar
Matthias Redies committed
113
114
115
      ! calculate radial basisfunctions u and u' with
      ! the spherical part of the potential vr0 and store them in
      ! bas1 = large component ,bas2 = small component
116

117
118
      call gen_bas_fun(atoms, enpara, gridf, input, hybrid, mpi, vr0, usdus, bas1, bas2)

Matthias Redies's avatar
Matthias Redies committed
119
120
121
      ! - - - - - - SETUP OF THE MIXED BASIS IN THE IR - - - - - - -

      ! construct G-vectors with cutoff smaller than gcutm
Matthias Redies's avatar
Matthias Redies committed
122
      call mpbasis%gen_gvec(cell, kpts, mpi)
Matthias Redies's avatar
Matthias Redies committed
123
124
125
126
127
128
129
130

      ! - - - - - - - - Set up MT product basis for the non-local exchange potential  - - - - - - - - - -

      IF (mpi%irank == 0) THEN
         WRITE (6, '(A)') 'MT product basis for non-local exchange potential:'
         WRITE (6, '(A)') 'Reduction due to overlap (quality of orthonormality, should be < 1.0E-06)'
      END IF

131
      allocate(mpbasis%num_radbasfn(0:maxval(hybrid%lcutm1), atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
132
133
      allocate(seleco(maxval(hybrid%num_radfun_per_l), 0:atoms%lmaxd))
      allocate(selecu(maxval(hybrid%num_radfun_per_l), 0:atoms%lmaxd))
134
      mpbasis%num_radbasfn = 0    !!! 01/12/10 jij%M.b.
Matthias Redies's avatar
Matthias Redies committed
135

136
      ! determine maximal indices of (radial) mixed-basis functions (->num_radbasfn)
Matthias Redies's avatar
Matthias Redies committed
137
      ! (will be reduced later-on due to overlap)
138
      hybrid%max_indx_p_1 = 0
Matthias Redies's avatar
Matthias Redies committed
139
140
141
142
143
144
145
146
147
      DO itype = 1, atoms%ntype
         seleco = .FALSE.
         selecu = .FALSE.
         seleco(1, 0:hybrid%select1(1, itype)) = .TRUE.
         selecu(1, 0:hybrid%select1(3, itype)) = .TRUE.
         seleco(2, 0:hybrid%select1(2, itype)) = .TRUE.
         selecu(2, 0:hybrid%select1(4, itype)) = .TRUE.

         ! include local orbitals
Matthias Redies's avatar
Matthias Redies committed
148
         IF (maxval(hybrid%num_radfun_per_l) >= 3) THEN
Matthias Redies's avatar
Matthias Redies committed
149
150
            seleco(3:,:) = .TRUE.
            selecu(3:,:) = .TRUE.
Matthias Redies's avatar
Matthias Redies committed
151
152
153
         END IF

         DO l = 0, hybrid%lcutm1(itype)
Matthias Redies's avatar
Matthias Redies committed
154
            n_radbasfn = 0
Matthias Redies's avatar
Matthias Redies committed
155
156
157
158
159
            M = 0

            !
            ! valence * valence
            !
Matthias Redies's avatar
Matthias Redies committed
160
161
162
163
164
165
            if(.not. allocated(selecmat)) then
               allocate(selecmat(maxval(hybrid%num_radfun_per_l), &
                                 0:atoms%lmaxd, &
                                 maxval(hybrid%num_radfun_per_l), &
                                 0:atoms%lmaxd))
            endif
Matthias Redies's avatar
Matthias Redies committed
166
            selecmat = calc_selecmat(atoms, hybrid, seleco, selecu)
Matthias Redies's avatar
Matthias Redies committed
167
168
169
170

            DO l1 = 0, atoms%lmax(itype)
               DO l2 = 0, atoms%lmax(itype)
                  IF (l >= ABS(l1 - l2) .AND. l <= l1 + l2) THEN
171
172
                     DO n1 = 1, hybrid%num_radfun_per_l(l1, itype)
                        DO n2 = 1, hybrid%num_radfun_per_l(l2, itype)
Matthias Redies's avatar
Matthias Redies committed
173
174
                           M = M + 1
                           IF (selecmat(n1, l1, n2, l2)) THEN
Matthias Redies's avatar
Matthias Redies committed
175
                              n_radbasfn = n_radbasfn + 1
Matthias Redies's avatar
Matthias Redies committed
176
177
178
179
                              selecmat(n2, l2, n1, l1) = .FALSE. ! prevent double counting of products (a*b = b*a)
                           END IF
                        END DO
                     END DO
180
181
182
                  END IF
               END DO
            END DO
Matthias Redies's avatar
Matthias Redies committed
183
            IF (n_radbasfn == 0 .AND. mpi%irank == 0) &
Matthias Redies's avatar
Matthias Redies committed
184
185
               WRITE (6, '(A)') 'mixedbasis: Warning!  No basis-function product of '//lchar(l)// &
               '-angular momentum defined.'
186
            hybrid%max_indx_p_1 = MAX(hybrid%max_indx_p_1, M)
187
            mpbasis%num_radbasfn(l, itype) = n_radbasfn*input%jspins
188
189
         END DO
      END DO
Matthias Redies's avatar
Matthias Redies committed
190

Matthias Redies's avatar
Matthias Redies committed
191
      allocate(mpbasis%radbasfn_mt(atoms%jmtd,&
192
                            maxval(mpbasis%num_radbasfn), &
193
194
                            0:maxval(hybrid%lcutm1), &
                            atoms%ntype), source=0.0)
Matthias Redies's avatar
Matthias Redies committed
195
196
197
198
199
200
201
202
203
204
205

      ! Define product bases and reduce them according to overlap

      DO itype = 1, atoms%ntype
         seleco = .FALSE.
         selecu = .FALSE.
         seleco(1, 0:hybrid%select1(1, itype)) = .TRUE.
         selecu(1, 0:hybrid%select1(3, itype)) = .TRUE.
         seleco(2, 0:hybrid%select1(2, itype)) = .TRUE.
         selecu(2, 0:hybrid%select1(4, itype)) = .TRUE.
         ! include lo's
Matthias Redies's avatar
Matthias Redies committed
206
         IF (maxval(hybrid%num_radfun_per_l) >= 3) THEN
Matthias Redies's avatar
Matthias Redies committed
207
208
            seleco(3:,:) = .TRUE.
            selecu(3:,:) = .TRUE.
Matthias Redies's avatar
Matthias Redies committed
209
         END IF
Matthias Redies's avatar
Matthias Redies committed
210

Matthias Redies's avatar
Matthias Redies committed
211
         n_grid_pt = atoms%jri(itype)
Matthias Redies's avatar
Matthias Redies committed
212
         DO l = 0, hybrid%lcutm1(itype)
213
            n_radbasfn = mpbasis%num_radbasfn(l, itype)
Matthias Redies's avatar
Matthias Redies committed
214
215
            ! allow for zero product-basis functions for
            ! current l-quantum number
Matthias Redies's avatar
Matthias Redies committed
216
            IF (n_radbasfn == 0) THEN
Matthias Redies's avatar
Matthias Redies committed
217
218
219
220
221
               IF (mpi%irank == 0) WRITE (6, '(6X,A,'':   0 ->   0'')') lchar(l)
               CYCLE
            END IF

            ! set up the overlap matrix
222
223
224
            allocate(eigv(n_radbasfn, n_radbasfn), source=0.0)
            allocate(work(3*n_radbasfn), source=0.0)
            allocate(eig(n_radbasfn), source=0.0)
225
            i_basfn = 0
Matthias Redies's avatar
Matthias Redies committed
226
227

            ! valence*valence
Matthias Redies's avatar
Matthias Redies committed
228
            selecmat =  calc_selecmat(atoms, hybrid, seleco, selecu)
Matthias Redies's avatar
Matthias Redies committed
229
230
231

            DO l1 = 0, atoms%lmax(itype)
               DO l2 = 0, atoms%lmax(itype)
Matthias Redies's avatar
Matthias Redies committed
232
233
234
                  IF (l >= ABS(l1 - l2) .AND. l <= l1 + l2) THEN
                     DO n1 = 1, hybrid%num_radfun_per_l(l1, itype)
                        DO n2 = 1, hybrid%num_radfun_per_l(l2, itype)
Matthias Redies's avatar
Matthias Redies committed
235

Matthias Redies's avatar
Matthias Redies committed
236
237
                           IF (selecmat(n1, l1, n2, l2)) THEN
                              DO jspin = 1, input%jspins
238
239
                                 i_basfn = i_basfn + 1
                                 IF (i_basfn > n_radbasfn) call judft_error('got too many product functions', hint='This is a BUG, please report', calledby='mixedbasis')
Matthias Redies's avatar
Matthias Redies committed
240

Matthias Redies's avatar
Matthias Redies committed
241
                                 mpbasis%radbasfn_mt(:n_grid_pt, i_basfn, l, itype) &
242
                                    = (   bas1(:n_grid_pt, n1, l1, itype, jspin) &
Matthias Redies's avatar
Matthias Redies committed
243
244
245
246
                                        * bas1(:n_grid_pt, n2, l2, itype, jspin) &
                                        + bas2(:n_grid_pt, n1, l1, itype, jspin) &
                                        * bas2(:n_grid_pt, n2, l2, itype, jspin) &
                                      ) / atoms%rmsh(:n_grid_pt, itype)
Matthias Redies's avatar
Matthias Redies committed
247

Matthias Redies's avatar
Matthias Redies committed
248
                                 !normalize radbasfn_mt
Matthias Redies's avatar
Matthias Redies committed
249
250
251
                                 norm = calc_radbas_norm(atoms, mpbasis, l, i_basfn, itype, gridf)
                                 mpbasis%radbasfn_mt(:n_grid_pt, i_basfn, l, itype) &
                                 = mpbasis%radbasfn_mt(:n_grid_pt, i_basfn, l, itype)/norm
Matthias Redies's avatar
Matthias Redies committed
252

Matthias Redies's avatar
Matthias Redies committed
253
254
255
256
257
258
259
                              END DO !jspin
                              ! prevent double counting of products (a*b = b*a)
                              selecmat(n2, l2, n1, l1) = .FALSE.
                           END IF
                        END DO !n2
                     END DO !n1
                  ENDIF
Matthias Redies's avatar
Matthias Redies committed
260
261
262
               END DO !l2
            END DO  !l1

263
            IF (i_basfn /= n_radbasfn) call judft_error('counting error for product functions', hint='This is a BUG, please report', calledby='mixedbasis')
Matthias Redies's avatar
Matthias Redies committed
264

Matthias Redies's avatar
typo    
Matthias Redies committed
265
            ! In order to get rid of the linear dependencies in the
Matthias Redies's avatar
Matthias Redies committed
266
            ! radial functions radbasfn_mt belonging to fixed l and itype
Matthias Redies's avatar
Matthias Redies committed
267
268
269
270
            ! the overlap matrix is diagonalized and those eigenvectors
            ! with a eigenvalue greater then hybrid%tolerance1 are retained

            ! Calculate overlap
Matthias Redies's avatar
Matthias Redies committed
271
            call calc_olap_radbasfn(atoms, mpbasis, l, itype, gridf, olap)
Matthias Redies's avatar
Matthias Redies committed
272
273
274
275

            ! Diagonalize
            CALL diagonalize(eigv, eig, olap)

276
            call filter_radbasfn(hybrid, l, itype, n_radbasfn, eig, eigv, mpbasis)
Matthias Redies's avatar
Matthias Redies committed
277

Matthias Redies's avatar
Matthias Redies committed
278
            call trafo_to_orthonorm_bas(mpbasis, n_radbasfn, n_grid_pt, l, itype, eig, eigv)
279
            nn = mpbasis%num_radbasfn(l, itype)
Matthias Redies's avatar
Matthias Redies committed
280
281

            ! Add constant function to l=0 basis and then do a Gram-Schmidt orthonormalization
Matthias Redies's avatar
Matthias Redies committed
282
            call add_l0_fun(atoms, hybrid, n_grid_pt, l, itype, gridf, nn, mpbasis)
Matthias Redies's avatar
Matthias Redies committed
283
284

            ! Check orthonormality of product basis
285
            call mpbasis%check_orthonormality(atoms, mpi, l, itype, gridf)
Matthias Redies's avatar
Matthias Redies committed
286

287
            deallocate(olap, eigv, work, eig)
Matthias Redies's avatar
Matthias Redies committed
288
289

         END DO !l
290
         IF (mpi%irank == 0) WRITE (6, '(6X,A,I7)') 'Total:', SUM(mpbasis%num_radbasfn(0:hybrid%lcutm1(itype), itype))
Matthias Redies's avatar
Matthias Redies committed
291
292
      END DO ! itype

293
294
295
      allocate(basmhlp(atoms%jmtd, maxval(mpbasis%num_radbasfn), 0:maxval(hybrid%lcutm1), atoms%ntype))
      basmhlp(1:atoms%jmtd, 1:maxval(mpbasis%num_radbasfn), 0:maxval(hybrid%lcutm1), 1:atoms%ntype) &
         = mpbasis%radbasfn_mt(1:atoms%jmtd, 1:maxval(mpbasis%num_radbasfn), 0:maxval(hybrid%lcutm1), 1:atoms%ntype)
Matthias Redies's avatar
Matthias Redies committed
296
      deallocate(mpbasis%radbasfn_mt)
297
      allocate(mpbasis%radbasfn_mt(atoms%jmtd, maxval(mpbasis%num_radbasfn), 0:maxval(hybrid%lcutm1), atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
298
      mpbasis%radbasfn_mt = basmhlp
Matthias Redies's avatar
Matthias Redies committed
299

300
      deallocate(basmhlp, seleco, selecu, selecmat)
Matthias Redies's avatar
Matthias Redies committed
301
302
303
304
305
306
307
308
309
310
311
312
313
314

      !
      ! now we build linear combinations of the radial functions
      ! such that they possess no moment except one radial function in each l-channel
      !
      IF (mpi%irank == 0) THEN
         WRITE (6, '(/,A,/,A)') 'Build linear combinations of radial '// &
            'functions in each l-channel,', &
            'such that they possess no multipolmoment'// &
            ' except the last function:'

         WRITE (6, '(/,17x,A)') 'moment  (quality of orthonormality)'
      END IF
      DO itype = 1, atoms%ntype
Matthias Redies's avatar
Matthias Redies committed
315
         n_grid_pt = atoms%jri(itype)
Matthias Redies's avatar
Matthias Redies committed
316
317
318
319
320
321

         IF (atoms%ntype > 1 .AND. mpi%irank == 0) WRITE (6, '(6X,A,I3)') 'Atom type', itype

         DO l = 0, hybrid%lcutm1(itype)
            ! determine radial function with the largest moment
            ! this function is used to build the linear combinations
322
            max_momentum = 0
323
            DO i = 1, mpbasis%num_radbasfn(l, itype)
324
               momentum = intgrf(atoms%rmsh(:n_grid_pt, itype)**(l + 1)*mpbasis%radbasfn_mt(:n_grid_pt, i, l, itype), &
Matthias Redies's avatar
Matthias Redies committed
325
                              atoms, itype, gridf)
326
               IF (ABS(momentum) > max_momentum) THEN
Matthias Redies's avatar
Matthias Redies committed
327
                  n_radbasfn = i
328
                  max_momentum = momentum
Matthias Redies's avatar
Matthias Redies committed
329
330
331
332
               END IF
            END DO

            ! rearrange order of radial functions such that the last function possesses the largest moment
Matthias Redies's avatar
Matthias Redies committed
333
            bashlp(:n_grid_pt) = mpbasis%radbasfn_mt(:n_grid_pt, n_radbasfn, l, itype)
Matthias Redies's avatar
Matthias Redies committed
334
            mpbasis%radbasfn_mt(:n_grid_pt,&
335
                                n_radbasfn:mpbasis%num_radbasfn(l, itype)-1,&
Matthias Redies's avatar
Matthias Redies committed
336
337
                                :, itype)&
               =  mpbasis%radbasfn_mt(:n_grid_pt,&
338
                                      n_radbasfn+1:mpbasis%num_radbasfn(l, itype),&
Matthias Redies's avatar
Matthias Redies committed
339
340
                                      :, itype)
            mpbasis%radbasfn_mt(:n_grid_pt, &
341
                                mpbasis%num_radbasfn(l, itype),&
Matthias Redies's avatar
Matthias Redies committed
342
343
                                l, itype) &
               = bashlp(:n_grid_pt)
Matthias Redies's avatar
Matthias Redies committed
344
345
346
347
348
         END DO

         DO l = 0, hybrid%lcutm1(itype)
            IF (mpi%irank == 0) WRITE (6, '(6X,A)') lchar(l)//':'

349
            IF (mpbasis%num_radbasfn(l, itype) == 0) THEN
Matthias Redies's avatar
Matthias Redies committed
350
351
352
353
               IF (mpi%irank == 0) WRITE (6, '(6X,A,'':   0 ->    '')') lchar(l)
               CYCLE
            END IF

354
            n_radbasfn = mpbasis%num_radbasfn(l, itype)
Matthias Redies's avatar
Matthias Redies committed
355
            DO i = 1, n_radbasfn-1
Matthias Redies's avatar
Matthias Redies committed
356
               ! calculate moment of radial function i
Matthias Redies's avatar
Matthias Redies committed
357
               rdum1 = intgrf(atoms%rmsh(:n_grid_pt, itype)**(l + 1)*mpbasis%radbasfn_mt(:n_grid_pt, i, l, itype), &
Matthias Redies's avatar
Matthias Redies committed
358
                              atoms, itype, gridf)
359

Matthias Redies's avatar
Matthias Redies committed
360
               rdum = intgrf(atoms%rmsh(:n_grid_pt, itype)**(l + 1)*mpbasis%radbasfn_mt(:n_grid_pt, n_radbasfn, l, itype), &
Matthias Redies's avatar
Matthias Redies committed
361
                             atoms, itype, gridf)
Matthias Redies's avatar
Matthias Redies committed
362

Matthias Redies's avatar
Matthias Redies committed
363
               bashlp(:n_grid_pt) = mpbasis%radbasfn_mt(:n_grid_pt, n_radbasfn, l, itype)
Matthias Redies's avatar
Matthias Redies committed
364
365

               IF (SQRT(rdum**2 + rdum1**2) <= 1E-06 .AND. mpi%irank == 0) &
Matthias Redies's avatar
typo    
Matthias Redies committed
366
                  WRITE (6, *) 'Warning: Norm is smaller than 1E-06!'
Matthias Redies's avatar
Matthias Redies committed
367

Matthias Redies's avatar
Matthias Redies committed
368
               ! change function n_radbasfn such that n_radbasfn is orthogonal to i
Matthias Redies's avatar
Matthias Redies committed
369
               ! since the functions radbasfn_mt have been orthogonal on input
Matthias Redies's avatar
Matthias Redies committed
370
               ! the linear combination does not destroy the orthogonality to the residual functions
Matthias Redies's avatar
Matthias Redies committed
371
372
               mpbasis%radbasfn_mt(:n_grid_pt, n_radbasfn, l, itype) = rdum/SQRT(rdum**2 + rdum1**2)*bashlp(:n_grid_pt) &
                                                + rdum1/SQRT(rdum**2 + rdum1**2)*mpbasis%radbasfn_mt(:n_grid_pt, i, l, itype)
Matthias Redies's avatar
Matthias Redies committed
373

Matthias Redies's avatar
Matthias Redies committed
374
               ! combine basis function i and n_radbasfn so that they possess no momemt
Matthias Redies's avatar
Matthias Redies committed
375
376
               mpbasis%radbasfn_mt(:n_grid_pt, i, l, itype) = rdum1/SQRT(rdum**2 + rdum1**2)*bashlp(:n_grid_pt) &
                                                - rdum/SQRT(rdum**2 + rdum1**2)*mpbasis%radbasfn_mt(:n_grid_pt, i, l, itype)
Matthias Redies's avatar
Matthias Redies committed
377

Matthias Redies's avatar
Matthias Redies committed
378
               rdum1 = intgrf(atoms%rmsh(:n_grid_pt, itype)**(l + 1)*mpbasis%radbasfn_mt(:n_grid_pt, i, l, itype), &
Matthias Redies's avatar
Matthias Redies committed
379
                              atoms, itype, gridf)
Matthias Redies's avatar
Matthias Redies committed
380
381
382
383
384

               IF (rdum1 > 1E-10) call judft_error('moment of radial function does not vanish', calledby='mixedbasis')

               IF (mpi%irank == 0) WRITE (6, '(6x,I4,'' ->  '',ES8.1)') i, rdum1
            END DO
385

Matthias Redies's avatar
Matthias Redies committed
386
            ! test orthogonality
Matthias Redies's avatar
Matthias Redies committed
387
            call mpbasis%check_orthonormality(atoms, mpi, l, itype, gridf)
Matthias Redies's avatar
Matthias Redies committed
388
         ENDDO
389
      END DO
Matthias Redies's avatar
Matthias Redies committed
390

Matthias Redies's avatar
Matthias Redies committed
391
      call mpbasis%check_radbasfn(atoms, hybrid)
Matthias Redies's avatar
Matthias Redies committed
392
393
394
395
396
397

      !count basis functions
      hybrid%nbasp = 0
      DO itype = 1, atoms%ntype
         DO i = 1, atoms%neq(itype)
            DO l = 0, hybrid%lcutm1(itype)
398
               hybrid%nbasp = hybrid%nbasp + (2*l+1) * mpbasis%num_radbasfn(l, itype)
Matthias Redies's avatar
Matthias Redies committed
399
400
            END DO
         END DO
401
      END DO
Matthias Redies's avatar
Matthias Redies committed
402
      hybrid%maxbasm1 = hybrid%nbasp + maxval(mpbasis%ngptm)
Matthias Redies's avatar
Matthias Redies committed
403
      hybrid%nbasm = hybrid%nbasp + mpbasis%ngptm
Daniel Wortmann's avatar
Daniel Wortmann committed
404

Matthias Redies's avatar
Matthias Redies committed
405
406
      hybrid%maxlmindx = 0
      do itype = 1,atoms%ntype
Matthias Redies's avatar
Matthias Redies committed
407
         hybrid%maxlmindx = max(hybrid%maxlmindx,
Matthias Redies's avatar
Matthias Redies committed
408
409
                                SUM([hybrid%num_radfun_per_l(l, itype)*(2*l + 1), l=0, atoms%lmax(itype)]))
      enddo
Matthias Redies's avatar
Matthias Redies committed
410
   END SUBROUTINE mixedbasis
Matthias Redies's avatar
Matthias Redies committed
411

412
   subroutine gen_bas_fun(atoms, enpara, gridf, input, hybrid, mpi, vr0, usdus, bas1, bas2)
Matthias Redies's avatar
Matthias Redies committed
413
414
      use m_judft
      use m_types
415
416
      USE m_radfun, ONLY: radfun
      USE m_radflo, ONLY: radflo
Matthias Redies's avatar
Matthias Redies committed
417
      USE m_intgrf,   ONLY: intgrf
Matthias Redies's avatar
Matthias Redies committed
418
      implicit NONE
419
420
421
422
423
424
425
426
      type(t_atoms), intent(in)        :: atoms
      type(t_enpara), intent(in)       :: enpara
      type(t_input), intent(in)        :: input
      type(t_hybrid), intent(in)       :: hybrid
      type(t_mpi), intent(in)          :: mpi
      type(t_usdus), intent(inout)     :: usdus

      REAL, ALLOCATABLE, INTENT(INOUT) :: bas1(:,:,:,:,:), bas2(:,:,:,:,:)
427
      REAL, intent(in)                 :: vr0(:,:,:), gridf(:,:)
428
429
430
431
432
433
434
435

      REAL    ::   u(atoms%jmtd, 2, 0:atoms%lmaxd)
      REAL    ::  du(atoms%jmtd, 2, 0:atoms%lmaxd)
      REAL    :: flo(atoms%jmtd, 2, atoms%nlod)

      REAL    :: uuilon(atoms%nlod, atoms%ntype)
      REAL    :: duilon(atoms%nlod, atoms%ntype)
      REAL    :: ulouilopn(atoms%nlod, atoms%nlod, atoms%ntype)
436
      REAL    :: wronk, norm
Matthias Redies's avatar
Matthias Redies committed
437

438
      INTEGER :: itype, jspin, i, l, ilo, ok
Matthias Redies's avatar
Matthias Redies committed
439
      INTEGER :: n_grid_pt, noded, nodem
Matthias Redies's avatar
Matthias Redies committed
440
441
442
443
      u  = 0.0
      du = 0.0

      ! this is 5-D array. it could cause Problems in bigger systems
444
      allocate(bas1(atoms%jmtd,    &
Matthias Redies's avatar
Matthias Redies committed
445
                    maxval(hybrid%num_radfun_per_l), &
446
447
                    0:atoms%lmaxd, &
                    atoms%ntype,   &
Matthias Redies's avatar
Matthias Redies committed
448
449
450
451
452
453
454
                    input%jspins),   source=0.0, stat=ok)
      if(ok /= 0) call judft_error("Can't allocate bas1 array. Stat= " // int2str(ok))

      allocate(bas2, source=bas1, stat=ok)
      if(ok /= 0) call judft_error("Can't allocate bas1 array. Stat= " // int2str(ok))

      DO itype = 1, atoms%ntype
Matthias Redies's avatar
Matthias Redies committed
455
         n_grid_pt = atoms%jri(itype) ! number of radial gridpoints
Matthias Redies's avatar
Matthias Redies committed
456
457
458
459
460
         DO jspin = 1, input%jspins
            DO l = 0, atoms%lmax(itype)
               CALL radfun(l, itype, jspin, enpara%el0(l, itype, jspin), vr0(:,itype, jspin), atoms, &
                           u(:,:,l), du(:,:,l), usdus, nodem, noded, wronk)
            END DO
Matthias Redies's avatar
Matthias Redies committed
461
462
463
464
            bas1(1:n_grid_pt, 1, 0:atoms%lmaxd, itype, jspin)  = u(1:n_grid_pt, 1, 0:atoms%lmaxd)
            bas2(1:n_grid_pt, 1, 0:atoms%lmaxd, itype, jspin)  = u(1:n_grid_pt, 2, 0:atoms%lmaxd)
            bas1(1:n_grid_pt, 2, 0:atoms%lmaxd, itype, jspin) = du(1:n_grid_pt, 1, 0:atoms%lmaxd)
            bas2(1:n_grid_pt, 2, 0:atoms%lmaxd, itype, jspin) = du(1:n_grid_pt, 2, 0:atoms%lmaxd)
Matthias Redies's avatar
Matthias Redies committed
465
466
467
468
469
470
471

            ! generate radial functions for local orbitals
            IF (atoms%nlo(itype) >= 1) THEN
               CALL radflo(atoms, itype, jspin, enpara%ello0(1, 1, jspin), vr0(:,itype, jspin), &
                           u, du, mpi, usdus, uuilon, duilon, ulouilopn, flo)

               DO ilo = 1, atoms%nlo(itype)
Matthias Redies's avatar
Matthias Redies committed
472
473
                  bas1(1:n_grid_pt, 2+ilo, atoms%llo(ilo, itype), itype, jspin) = flo(1:n_grid_pt, 1, ilo)
                  bas2(1:n_grid_pt, 2+ilo, atoms%llo(ilo, itype), itype, jspin) = flo(1:n_grid_pt, 2, ilo)
Matthias Redies's avatar
Matthias Redies committed
474
475
476
477
               END DO
            END IF
         END DO
      END DO
478
479
480
481
482
483
484

      ! the radial functions are normalized
      DO jspin = 1, input%jspins
         DO itype = 1, atoms%ntype
            DO l = 0, atoms%lmax(itype)
               DO i = 1, hybrid%num_radfun_per_l(l, itype)
                  norm = sqrt(intgrf(bas1(:,i, l, itype, jspin)**2 + bas2(:,i, l, itype, jspin)**2, &
Matthias Redies's avatar
Matthias Redies committed
485
                                atoms, itype, gridf))
486
487
488
489
490
491
                  bas1(:atoms%jri(itype), i, l, itype, jspin) = bas1(:atoms%jri(itype), i, l, itype, jspin)/norm
                  bas2(:atoms%jri(itype), i, l, itype, jspin) = bas2(:atoms%jri(itype), i, l, itype, jspin)/norm
               END DO
            END DO
         END DO
      END DO
Matthias Redies's avatar
Matthias Redies committed
492
   end subroutine gen_bas_fun
493

Matthias Redies's avatar
Matthias Redies committed
494
   function calc_selecmat(atoms,hybrid,seleco, selecu) result(selecmat)
Matthias Redies's avatar
Matthias Redies committed
495
      ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
Matthias Redies's avatar
Matthias Redies committed
496
      use m_types
Matthias Redies's avatar
Matthias Redies committed
497
      use m_judft
498
      implicit NONE
Matthias Redies's avatar
Matthias Redies committed
499

Matthias Redies's avatar
Matthias Redies committed
500
501
502
503
504
505
      type(t_atoms),  intent(in) :: atoms
      type(t_hybrid), intent(in) :: hybrid
      LOGICAL, intent(in) :: seleco(maxval(hybrid%num_radfun_per_l), 0:atoms%lmaxd)
      LOGICAL, intent(in) :: selecu(maxval(hybrid%num_radfun_per_l), 0:atoms%lmaxd)
      LOGICAL  ::  selecmat(maxval(hybrid%num_radfun_per_l), 0:atoms%lmaxd, &
                            maxval(hybrid%num_radfun_per_l), 0:atoms%lmaxd)
506
507
508
509
      integer                       :: n1, l1, n2, l2

      ! column-major means left-most index varies the fastest
      do l2=0,atoms%lmaxd
Matthias Redies's avatar
Matthias Redies committed
510
         do n2=1,maxval(hybrid%num_radfun_per_l)
511
            do l1=0,atoms%lmaxd
Matthias Redies's avatar
Matthias Redies committed
512
               do n1=1,maxval(hybrid%num_radfun_per_l)
513
514
515
516
517
518
                  selecmat(n1,l1,n2,l2) = seleco(n1, l1) .AND. selecu(n2, l2)
               enddo
            enddo
         enddo
      enddo
   end function calc_selecmat
Matthias Redies's avatar
Matthias Redies committed
519

Matthias Redies's avatar
Matthias Redies committed
520
   function calc_radbas_norm(atoms, mpbasis, l, i_basfn, itype, gridf) result(norm)
521
      USE m_intgrf, ONLY: intgrf
Matthias Redies's avatar
Matthias Redies committed
522
523
      use m_types
      implicit NONE
Matthias Redies's avatar
Matthias Redies committed
524
525
526
527
      type(t_atoms), intent(in)   :: atoms
      type(t_mpbasis), intent(in) :: mpbasis
      integer, intent(in)         :: l, i_basfn, itype
      real, intent(in)            :: gridf(:,:)
Matthias Redies's avatar
Matthias Redies committed
528

Matthias Redies's avatar
Matthias Redies committed
529
      real                        :: norm
Matthias Redies's avatar
Matthias Redies committed
530

531
      norm = SQRT( &
Matthias Redies's avatar
Matthias Redies committed
532
                  intgrf(mpbasis%radbasfn_mt(:,i_basfn, l, itype)**2, &
Matthias Redies's avatar
Matthias Redies committed
533
                         atoms,&
Matthias Redies's avatar
Matthias Redies committed
534
535
536
                         itype, gridf)&
                 )
   end function calc_radbas_norm
Matthias Redies's avatar
Matthias Redies committed
537

538
539
540
541
542
543
544

   subroutine filter_radbasfn(hybrid, l, itype, n_radbasfn, eig, eigv, mpbasis)

                  ! Get rid of linear dependencies (eigenvalue <= hybrid%tolerance1)
      use m_types
      implicit none
      type(t_hybrid), intent(in)     :: hybrid
545
546
      integer, intent(in)          :: l, itype, n_radbasfn
      real, intent(inout)          :: eig(:), eigv(:,:)
547
548
549
550
551
552
553
554
      type(t_mpbasis), intent(inout) :: mpbasis

      integer              :: num_radbasfn, i_bas
      integer, allocatable :: remaining_basfn(:)

      allocate(remaining_basfn(n_radbasfn), source=1)
      num_radbasfn = 0

555
      DO i_bas = 1, mpbasis%num_radbasfn(l, itype)
556
557
558
559
560
561
         IF (eig(i_bas) > hybrid%tolerance1) THEN
            num_radbasfn = num_radbasfn + 1
            remaining_basfn(num_radbasfn) = i_bas
         END IF
      END DO

562
      mpbasis%num_radbasfn(l, itype) = num_radbasfn
563
564
565
      eig = eig(remaining_basfn)
      eigv(:,:) = eigv(:,remaining_basfn)
   end subroutine filter_radbasfn
566

Matthias Redies's avatar
Matthias Redies committed
567
   subroutine trafo_to_orthonorm_bas(mpbasis, n_radbasfn, n_grid_pt, l, itype, eig, eigv)
568
569
      use m_types
      implicit NONE
Matthias Redies's avatar
Matthias Redies committed
570
      type(t_mpbasis), intent(inout)   :: mpbasis
571
572
573
574
575
      integer, intent(in)           :: n_radbasfn, n_grid_pt, l, itype
      real, intent(in)              :: eig(:), eigv(:,:)

      integer :: nn, i

576
      nn = mpbasis%num_radbasfn(l, itype)
577
      DO i = 1, n_grid_pt
Matthias Redies's avatar
Matthias Redies committed
578
579
         mpbasis%radbasfn_mt(i, 1:nn, l, itype) &
            = MATMUL(mpbasis%radbasfn_mt(i, 1:n_radbasfn, l, itype), eigv(:,1:nn))/SQRT(eig(:nn))
580
581
      END DO
   end subroutine trafo_to_orthonorm_bas
Matthias Redies's avatar
Matthias Redies committed
582
583
584
585
586
587
588
589
590
591
592
593
594
595

   subroutine add_l0_fun(atoms, hybrid, n_grid_pt, l, itype, gridf, nn, mpbasis)
      use m_types
      USE m_intgrf, ONLY: intgrf
      implicit none
      type(t_atoms), intent(in)       :: atoms
      type(t_hybrid), intent(in)      :: hybrid
      integer, intent(in)             :: n_grid_pt, l, itype
      real, intent(in)                :: gridf(:,:)
      integer, intent(inout)          :: nn
      type(t_mpbasis), intent(inout)  :: mpbasis

      integer                         :: i, j
      REAL, ALLOCATABLE               :: basmhlp(:,:,:,:)
596
      real                            :: norm
Matthias Redies's avatar
Matthias Redies committed
597
598
599
600
601
602
603
604
605
606
607
608
609

      IF (l == 0) THEN

         ! Check if radbasfn_mt must be reallocated
         IF (nn + 1 > SIZE(mpbasis%radbasfn_mt, 2)) THEN
            allocate(basmhlp(atoms%jmtd, nn + 1, 0:maxval(hybrid%lcutm1), atoms%ntype))
            basmhlp(:,1:nn, :,:) = mpbasis%radbasfn_mt
            deallocate(mpbasis%radbasfn_mt)
            allocate(mpbasis%radbasfn_mt(atoms%jmtd, nn + 1, 0:maxval(hybrid%lcutm1), atoms%ntype))
            mpbasis%radbasfn_mt(:,1:nn, :,:) = basmhlp(:,1:nn, :,:)
            deallocate(basmhlp)
         END IF

610
611
612
613
614
         ! add l = 0 function
         mpbasis%radbasfn_mt(:n_grid_pt, nn + 1, 0, itype) &
           = atoms%rmsh(:n_grid_pt, itype) / SQRT(atoms%rmsh(n_grid_pt, itype)**3/3)

         ! perform gram-schmidt orthonormalization
Matthias Redies's avatar
Matthias Redies committed
615
616
         DO i = nn, 1, -1
            DO j = i + 1, nn + 1
617
618
619
620
621
622
623
               mpbasis%radbasfn_mt(:n_grid_pt, i, 0, itype) &
                  = mpbasis%radbasfn_mt(:n_grid_pt, i, 0, itype) &
                     - intgrf( &
                                mpbasis%radbasfn_mt(:n_grid_pt, i, 0, itype) &
                              * mpbasis%radbasfn_mt(:n_grid_pt, j, 0, itype), &
                              atoms, itype, gridf) &
                       * mpbasis%radbasfn_mt(:n_grid_pt, j, 0, itype)
Matthias Redies's avatar
Matthias Redies committed
624
625

            END DO
626
627
628
629
630

            ! renormalize
            norm = SQRT(intgrf(mpbasis%radbasfn_mt(:n_grid_pt, i, 0, itype)**2, atoms, itype, gridf))
            mpbasis%radbasfn_mt(:n_grid_pt, i, 0, itype) &
               = mpbasis%radbasfn_mt(:n_grid_pt, i, 0, itype) / norm
Matthias Redies's avatar
Matthias Redies committed
631
632
         END DO
         nn = nn + 1
633
         mpbasis%num_radbasfn(l, itype) = nn
Matthias Redies's avatar
Matthias Redies committed
634
635
      END IF
   end subroutine add_l0_fun
636
637


Daniel Wortmann's avatar
Daniel Wortmann committed
638
END MODULE m_mixedbasis