mixedbasis.F90 26.9 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
            ! the overlap matrix is diagonalized and those eigenvectors
268
            ! with a eigenvalue greater then mpbasis%linear_dep_tol are retained
Matthias Redies's avatar
Matthias Redies committed
269
270

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

            ! Diagonalize
274
            call mpbasis_diagonialize_olap(olap, eig, eigv)
Matthias Redies's avatar
Matthias Redies committed
275

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

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
408
409
         hybrid%maxlmindx = max(hybrid%maxlmindx,&
                                SUM([(hybrid%num_radfun_per_l(l, itype)*(2*l + 1), l=0, atoms%lmax(itype))])&
                                )
Matthias Redies's avatar
Matthias Redies committed
410
      enddo
Matthias Redies's avatar
Matthias Redies committed
411
   END SUBROUTINE mixedbasis
Matthias Redies's avatar
Matthias Redies committed
412

413
   subroutine gen_bas_fun(atoms, enpara, gridf, input, hybrid, mpi, vr0, usdus, bas1, bas2)
Matthias Redies's avatar
Matthias Redies committed
414
415
      use m_judft
      use m_types
416
417
      USE m_radfun, ONLY: radfun
      USE m_radflo, ONLY: radflo
Matthias Redies's avatar
Matthias Redies committed
418
      USE m_intgrf,   ONLY: intgrf
Matthias Redies's avatar
Matthias Redies committed
419
      implicit NONE
420
421
422
423
424
425
426
427
      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(:,:,:,:,:)
428
      REAL, intent(in)                 :: vr0(:,:,:), gridf(:,:)
429
430
431
432
433
434
435
436

      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)
437
      REAL    :: wronk, norm
Matthias Redies's avatar
Matthias Redies committed
438

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

      ! this is 5-D array. it could cause Problems in bigger systems
445
      allocate(bas1(atoms%jmtd,    &
Matthias Redies's avatar
Matthias Redies committed
446
                    maxval(hybrid%num_radfun_per_l), &
447
448
                    0:atoms%lmaxd, &
                    atoms%ntype,   &
Matthias Redies's avatar
Matthias Redies committed
449
450
451
452
453
454
455
                    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
456
         n_grid_pt = atoms%jri(itype) ! number of radial gridpoints
Matthias Redies's avatar
Matthias Redies committed
457
458
459
460
461
         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
462
463
464
465
            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
466
467
468
469
470
471
472

            ! 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
473
474
                  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
475
476
477
478
               END DO
            END IF
         END DO
      END DO
479
480
481
482
483
484
485

      ! 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
486
                                atoms, itype, gridf))
487
488
489
490
491
492
                  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
493
   end subroutine gen_bas_fun
494

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

Matthias Redies's avatar
Matthias Redies committed
501
502
503
504
505
506
      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)
507
508
509
510
      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
511
         do n2=1,maxval(hybrid%num_radfun_per_l)
512
            do l1=0,atoms%lmaxd
Matthias Redies's avatar
Matthias Redies committed
513
               do n1=1,maxval(hybrid%num_radfun_per_l)
514
515
516
517
518
519
                  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
520

Matthias Redies's avatar
Matthias Redies committed
521
   function calc_radbas_norm(atoms, mpbasis, l, i_basfn, itype, gridf) result(norm)
522
      USE m_intgrf, ONLY: intgrf
Matthias Redies's avatar
Matthias Redies committed
523
524
      use m_types
      implicit NONE
Matthias Redies's avatar
Matthias Redies committed
525
526
527
528
      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
529

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

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

539

540
   subroutine trafo_to_orthonorm_bas(mpbasis, n_radbasfn, n_grid_pt, l, itype, eig, eigv)
541
542
      use m_types
      implicit NONE
Matthias Redies's avatar
Matthias Redies committed
543
      type(t_mpbasis), intent(inout)   :: mpbasis
544
545
      integer, intent(in)           :: n_radbasfn, n_grid_pt, l, itype
      real, intent(in)              :: eig(:), eigv(:,:)
546
547
548

      integer :: nn, i

549
      nn = mpbasis%num_radbasfn(l, itype)
550
      DO i = 1, n_grid_pt
Matthias Redies's avatar
Matthias Redies committed
551
         mpbasis%radbasfn_mt(i, 1:nn, l, itype) &
552
            = MATMUL(mpbasis%radbasfn_mt(i, 1:n_radbasfn, l, itype), eigv(:,1:nn))/SQRT(eig(:nn))
553
554
      END DO
   end subroutine trafo_to_orthonorm_bas
Matthias Redies's avatar
Matthias Redies committed
555
556
557
558
559
560
561
562
563
564
565
566
567
568

   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(:,:,:,:)
569
      real                            :: norm
Matthias Redies's avatar
Matthias Redies committed
570
571
572
573
574
575
576
577
578
579
580
581
582

      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

583
584
585
586
587
         ! 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
588
589
         DO i = nn, 1, -1
            DO j = i + 1, nn + 1
590
591
592
593
594
595
596
               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
597
598

            END DO
599
600
601
602
603

            ! 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
604
605
         END DO
         nn = nn + 1
606
         mpbasis%num_radbasfn(l, itype) = nn
Matthias Redies's avatar
Matthias Redies committed
607
608
      END IF
   end subroutine add_l0_fun
609
610


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