types_mpbasis.f90 16.3 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1
2
3
4
module m_types_mpbasis
   implicit none

   type t_mpbasis
5
6
7
8
      INTEGER, ALLOCATABLE   :: gptm(:,:) ! (3, num_gpts)
      INTEGER, ALLOCATABLE   :: ngptm(:)
      INTEGER, ALLOCATABLE   :: gptm_ptr(:,:)
      REAL                   :: g_cutoff
9
      INTEGER, ALLOCATABLE   :: num_radbasfn(:,:)
Matthias Redies's avatar
Matthias Redies committed
10
      REAL, ALLOCATABLE      :: radbasfn_mt(:,:,:,:)
11
      REAL                   :: linear_dep_tol  !only read in
Matthias Redies's avatar
Matthias Redies committed
12
   CONTAINS
Matthias Redies's avatar
Matthias Redies committed
13
14
15
16
17
18
19
20
      procedure :: num_gpts               => mpbasis_num_gpts
      procedure :: gen_gvec               => mpbasis_gen_gvec
      procedure :: check_orthonormality   => mpbasis_check_orthonormality
      procedure :: check_radbasfn         => mpbasis_check_radbasfn
      procedure :: calc_olap_radbasfn     => mpbasis_calc_olap_radbasfn
      procedure :: filter_radbasfn        => mpbasis_filter_radbasfn
      procedure :: trafo_to_orthonorm_bas => mpbasis_trafo_to_orthonorm_bas
      procedure :: add_l0_fun             => mpbasis_add_l0_fun
Matthias Redies's avatar
Matthias Redies committed
21
      procedure :: reduce_linear_dep      => mpbasis_reduce_linear_dep
Matthias Redies's avatar
Matthias Redies committed
22
   end type t_mpbasis
Matthias Redies's avatar
Matthias Redies committed
23
24
25
26
27
28
29
30
contains
   function mpbasis_num_gpts(mpbasis)
      implicit NONE
      class(t_mpbasis), intent(in) :: mpbasis
      integer    :: mpbasis_num_gpts

      mpbasis_num_gpts = size(mpbasis%gptm,dim=2)
   end function mpbasis_num_gpts
Matthias Redies's avatar
Matthias Redies committed
31

Matthias Redies's avatar
Matthias Redies committed
32
   subroutine mpbasis_gen_gvec(mpbasis, cell, kpts, mpi)
Matthias Redies's avatar
Matthias Redies committed
33
34
      use m_types_setup
      use m_types_kpts
Matthias Redies's avatar
Matthias Redies committed
35
      use m_types_mpi
Matthias Redies's avatar
Matthias Redies committed
36
37
38
      USE m_intgrf, ONLY: intgrf_init, intgrf
      use m_rorder, only: rorderpf
      implicit NONE
39
      class(t_mpbasis), intent(inout) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
40
41
      type(t_cell), intent(in)       :: cell
      type(t_kpts), intent(in)       :: kpts
Matthias Redies's avatar
Matthias Redies committed
42
      type(t_mpi), intent(in)        :: mpi
Matthias Redies's avatar
Matthias Redies committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64


      integer :: i, n, n1, n2, divconq
      integer :: x, y, z, ikpt, igpt
      integer :: g(3)
      real    :: longest_k, kvec(3)
      logical :: l_found_new_gpt, l_found_kg_in_sphere

      INTEGER, ALLOCATABLE            ::  unsrt_pgptm(:,:) ! unsorted pointers to g vectors
      REAL, ALLOCATABLE               ::  length_kg(:,:) ! length of the vectors k + G
      INTEGER, ALLOCATABLE            ::  ptr(:)

      allocate(mpbasis%ngptm(kpts%nkptf))

      mpbasis%ngptm = 0
      i = 0
      n = -1

      longest_k = MAXVAL([(norm2(MATMUL(kpts%bkf(:,ikpt), cell%bmat)), ikpt=1, kpts%nkptf)])

      ! a first run for the determination of the dimensions of the fields gptm,pgptm

Matthias Redies's avatar
Matthias Redies committed
65
      do
Matthias Redies's avatar
Matthias Redies committed
66
67
         n = n + 1
         l_found_new_gpt = .FALSE.
Matthias Redies's avatar
Matthias Redies committed
68
         do x = -n, n
Matthias Redies's avatar
Matthias Redies committed
69
            n1 = n - ABS(x)
Matthias Redies's avatar
Matthias Redies committed
70
            do y = -n1, n1
Matthias Redies's avatar
Matthias Redies committed
71
               n2 = n1 - ABS(y)
Matthias Redies's avatar
Matthias Redies committed
72
               do z = -n2, n2, MAX(2*n2, 1)
Matthias Redies's avatar
Matthias Redies committed
73
74
75
                  g = [x, y, z]
                  IF ((norm2(MATMUL(g, cell%bmat)) - longest_k) > mpbasis%g_cutoff) CYCLE
                  l_found_kg_in_sphere = .FALSE.
Matthias Redies's avatar
Matthias Redies committed
76
                  do ikpt = 1, kpts%nkptf
Matthias Redies's avatar
Matthias Redies committed
77
78
79
80
81
82
83
84
85
                     IF (norm2(MATMUL(kpts%bkf(:,ikpt) + g, cell%bmat)) <= mpbasis%g_cutoff) THEN
                        IF (.NOT. l_found_kg_in_sphere) THEN
                           i = i + 1
                           l_found_kg_in_sphere = .TRUE.
                        END IF

                        mpbasis%ngptm(ikpt) = mpbasis%ngptm(ikpt) + 1
                        l_found_new_gpt = .TRUE.
                     END IF
Matthias Redies's avatar
Matthias Redies committed
86
87
88
89
                  enddo ! k-loop
               enddo
            enddo
         enddo
Matthias Redies's avatar
Matthias Redies committed
90
         IF (.NOT. l_found_new_gpt) EXIT
Matthias Redies's avatar
Matthias Redies committed
91
      enddo
Matthias Redies's avatar
Matthias Redies committed
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

      allocate(mpbasis%gptm(3,i)) ! i = gptmd
      allocate(mpbasis%gptm_ptr(maxval(mpbasis%ngptm), kpts%nkptf))

      ! Allocate and initialize arrays needed for G vector ordering
      allocate(unsrt_pgptm(maxval(mpbasis%ngptm), kpts%nkptf))
      allocate(length_kG(maxval(mpbasis%ngptm), kpts%nkptf))

      mpbasis%gptm = 0
      mpbasis%gptm_ptr = 0
      mpbasis%ngptm = 0

      i = 0
      n = -1

      length_kG = 0
      unsrt_pgptm = 0

Matthias Redies's avatar
Matthias Redies committed
110
      do
Matthias Redies's avatar
Matthias Redies committed
111
112
         n = n + 1
         l_found_new_gpt = .FALSE.
Matthias Redies's avatar
Matthias Redies committed
113
         do x = -n, n
Matthias Redies's avatar
Matthias Redies committed
114
            n1 = n - ABS(x)
Matthias Redies's avatar
Matthias Redies committed
115
            do y = -n1, n1
Matthias Redies's avatar
Matthias Redies committed
116
               n2 = n1 - ABS(y)
Matthias Redies's avatar
Matthias Redies committed
117
               do z = -n2, n2, MAX(2*n2, 1)
Matthias Redies's avatar
Matthias Redies committed
118
119
120
                  g = [x, y, z]
                  IF ((norm2(MATMUL(g, cell%bmat)) - longest_k) > mpbasis%g_cutoff) CYCLE
                  l_found_kg_in_sphere = .FALSE.
Matthias Redies's avatar
Matthias Redies committed
121
                  do ikpt = 1, kpts%nkptf
Matthias Redies's avatar
Matthias Redies committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
                     kvec = kpts%bkf(:,ikpt)

                     IF (norm2(MATMUL(kvec + g, cell%bmat)) <= mpbasis%g_cutoff) THEN
                        IF (.NOT. l_found_kg_in_sphere) THEN
                           i = i + 1
                           mpbasis%gptm(:,i) = g
                           l_found_kg_in_sphere = .TRUE.
                        END IF

                        mpbasis%ngptm(ikpt) = mpbasis%ngptm(ikpt) + 1
                        l_found_new_gpt = .TRUE.

                        ! Position of the vector is saved as pointer
                        unsrt_pgptm(mpbasis%ngptm(ikpt), ikpt) = i
                        ! Save length of vector k + G for array sorting
                        length_kG(mpbasis%ngptm(ikpt), ikpt) = norm2(MATMUL(kvec + g, cell%bmat))
                     END IF
Matthias Redies's avatar
Matthias Redies committed
139
140
141
142
                  enddo
               enddo
            enddo
         enddo
Matthias Redies's avatar
Matthias Redies committed
143
         IF (.NOT. l_found_new_gpt) EXIT
Matthias Redies's avatar
Matthias Redies committed
144
      enddo
Matthias Redies's avatar
Matthias Redies committed
145
146

      ! Sort pointers in array, so that shortest |k+G| comes first
Matthias Redies's avatar
Matthias Redies committed
147
      do ikpt = 1, kpts%nkptf
Matthias Redies's avatar
Matthias Redies committed
148
149
150
151
152
153
         allocate(ptr(mpbasis%ngptm(ikpt)))
         ! Divide and conquer algorithm for arrays > 1000 entries
         divconq = MAX(0, INT(1.443*LOG(0.001*mpbasis%ngptm(ikpt))))
         ! create pointers which correspond to a sorted array
         CALL rorderpf(ptr, length_kG(1:mpbasis%ngptm(ikpt), ikpt), mpbasis%ngptm(ikpt), divconq)
         ! rearrange old pointers
Matthias Redies's avatar
Matthias Redies committed
154
         do igpt = 1, mpbasis%ngptm(ikpt)
Matthias Redies's avatar
Matthias Redies committed
155
            mpbasis%gptm_ptr(igpt, ikpt) = unsrt_pgptm(ptr(igpt), ikpt)
Matthias Redies's avatar
Matthias Redies committed
156
         enddo
Matthias Redies's avatar
Matthias Redies committed
157
         deallocate(ptr)
Matthias Redies's avatar
Matthias Redies committed
158
      enddo
Matthias Redies's avatar
Matthias Redies committed
159
160
161
162
163
164
165
166

      IF (mpi%irank == 0) THEN
         WRITE (6, '(/A)') 'Mixed basis'
         WRITE (6, '(A,I5)') 'Number of unique G-vectors: ', mpbasis%num_gpts()
         WRITE (6, *)
         WRITE (6, '(3x,A)') 'IR Plane-wave basis with cutoff of gcutm (mpbasis%g_cutoff/2*input%rkmax):'
         WRITE (6, '(5x,A,I5)') 'Maximal number of G-vectors:', maxval(mpbasis%ngptm)
      END IF
167
   end subroutine mpbasis_gen_gvec
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188

   subroutine mpbasis_check_orthonormality(mpbasis, atoms, mpi, l, itype, gridf)
      USE m_intgrf, ONLY: intgrf
      use m_types_setup
      use m_types_mpi
      use m_judft
      implicit none
      class(t_mpbasis)          :: mpbasis
      type(t_atoms), intent(in) :: atoms
      type(t_mpi), intent(in)   :: mpi
      integer, intent(in)       :: itype, l
      real, intent(in)          :: gridf(:,:)

      real              :: overlap, error, cum_err_sq
      integer           :: i, j, n_radbasfn

      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']


189
      n_radbasfn = mpbasis%num_radbasfn(l, itype)
190
191

      cum_err_sq = 0
Matthias Redies's avatar
Matthias Redies committed
192
193
      do i = 1, n_radbasfn
         do j = 1, i
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
            overlap = intgrf(mpbasis%radbasfn_mt(:,i, l, itype)*mpbasis%radbasfn_mt(:,j, l, itype), &
                           atoms, itype, gridf)

            IF (i == j) THEN
               error = ABS(1 - overlap)
               cum_err_sq = cum_err_sq + overlap**2
            ELSE
               error = ABS(overlap)
               cum_err_sq = cum_err_sq + 2*error**2
            END IF

            IF (error > 1e-6) THEN
               IF (mpi%irank == 0) THEN
                  WRITE (6, '(A)') 'mixedbasis: Bad orthonormality of ' &
                     //lchar(l)//'-product basis. Increase tolerance.'
                  WRITE (6, '(12X,A,F9.6,A,2(I3.3,A))') 'Deviation of', &
                     error, ' occurred for (', i, ',', j, ')-overlap.'
               END IF
               CALL judft_error("Bad orthonormality of product basis", &
                                 hint='Increase tolerance', &
                                 calledby='mixedbasis%check_orthonormality')
            END IF

Matthias Redies's avatar
Matthias Redies committed
217
218
         enddo
      enddo
219
220
221
222
223
224

      IF (mpi%irank == 0) THEN
         WRITE (6, '(6X,A,I4,''   ('',ES8.1,'' )'')') &
                    lchar(l) // ':', n_radbasfn, SQRT(cum_err_sq)/n_radbasfn
      END IF
   end subroutine mpbasis_check_orthonormality
Matthias Redies's avatar
Matthias Redies committed
225
226
227
228
229
230
231
232
233
234
235
236
237

   subroutine mpbasis_check_radbasfn(mpbasis, atoms, hybrid)
      use m_judft
      use m_types_hybrid
      use m_types_setup
      implicit none
      class(t_mpbasis), intent(in) :: mpbasis
      type(t_atoms), intent(in)    :: atoms
      type(t_hybrid), intent(in)   :: hybrid

      integer :: itype

      do itype = 1, atoms%ntype
238
239
         IF (ANY(mpbasis%num_radbasfn(0:hybrid%lcutm1(itype), itype) == 0)) THEN
            call judft_error('any mpbasis%num_radbasfn eq 0', calledby='mixedbasis')
Matthias Redies's avatar
Matthias Redies committed
240
241
242
         endif
      enddo
   end subroutine mpbasis_check_radbasfn
Matthias Redies's avatar
Matthias Redies committed
243
244
245

   subroutine mpbasis_calc_olap_radbasfn(mpbasis, atoms, l, itype, gridf, olap)
      USE m_intgrf, ONLY: intgrf
Matthias Redies's avatar
Matthias Redies committed
246
      use m_types_setup
Matthias Redies's avatar
Matthias Redies committed
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
      implicit NONE
      class(t_mpbasis), intent(in)       :: mpbasis
      type(t_atoms), intent(in)          :: atoms
      integer, intent(in)                :: l, itype
      real, intent(in)                   :: gridf(:,:)
      real, intent(inout), allocatable   :: olap(:,:)

      integer  :: n1, n2, n_radbasfn

      n_radbasfn = mpbasis%num_radbasfn(l, itype)
      if(allocated(olap)) then
         if(any(shape(olap) /= n_radbasfn)) then
            deallocate(olap)
         endif
      endif
      if(.not. allocated(olap)) allocate(olap(n_radbasfn, n_radbasfn), source=0.0)

      DO n2 = 1, n_radbasfn
         DO n1 = 1, n2
            olap(n1, n2) = intgrf(mpbasis%radbasfn_mt(:,n1, l, itype)*mpbasis%radbasfn_mt(:,n2, l, itype), &
                                  atoms, itype, gridf)
            olap(n2, n1) = olap(n1, n2)
         END DO
      END DO
   end subroutine mpbasis_calc_olap_radbasfn
Matthias Redies's avatar
Matthias Redies committed
272
273
274
275

   subroutine mpbasis_filter_radbasfn(mpbasis, l, itype, n_radbasfn, eig, eigv)
      ! Get rid of linear dependencies (eigenvalue <= mpbasis%linear_dep_tol)
      implicit none
276
277
278
      class(t_mpbasis), intent(inout)       :: mpbasis
      integer, intent(in)                   :: l, itype, n_radbasfn
      real, intent(inout)                   :: eig(:), eigv(:,:)
Matthias Redies's avatar
Matthias Redies committed
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296

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

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

      DO i_bas = 1, mpbasis%num_radbasfn(l, itype)
         IF (eig(i_bas) > mpbasis%linear_dep_tol) THEN
            num_radbasfn = num_radbasfn + 1
            remaining_basfn(num_radbasfn) = i_bas
         END IF
      END DO

      mpbasis%num_radbasfn(l, itype) = num_radbasfn
      eig = eig(remaining_basfn)
      eigv(:,:) = eigv(:,remaining_basfn)
   end subroutine mpbasis_filter_radbasfn
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314

   subroutine mpbasis_diagonialize_olap(olap, eig_val, eig_vec)
      use m_judft
      implicit NONE
      real, intent(in)  :: olap(:,:)
      real, allocatable :: eig_val(:), eig_vec(:,:)

      integer              :: n, size_iwork, info
      real                 :: size_work
      integer, allocatable :: iwork(:)
      real, allocatable    :: work(:)

      if(size(olap, dim=1) /= size(olap, dim=2)) then
         call juDFT_error("only square matrices can be diagonalized")
      endif

      n = size(olap, dim=1)

Matthias Redies's avatar
Matthias Redies committed
315
316
317
      if(allocated(eig_val)) then
         if(size(eig_val) /= n) deallocate(eig_val)
      endif
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
      if(.not. allocated(eig_val)) allocate(eig_val(n))

      eig_vec = olap
      ! get sizes of work arrays
      call dsyevd('V', 'U', n, eig_vec, n, eig_val,&
                  size_work, -1, size_iwork, -1, info)
      if(info /= 0) call juDFT_error("diagonalization for size failed")

      allocate(work(int(size_work)))
      allocate(iwork(size_iwork))

      call dsyevd('V', 'U', n, eig_vec, n, eig_val,&
                  work, int(size_work), iwork, size_iwork, info)
      if(info /= 0) call juDFT_error("diagonalization failed")

   end subroutine mpbasis_diagonialize_olap
Matthias Redies's avatar
Matthias Redies committed
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405

   subroutine mpbasis_trafo_to_orthonorm_bas(mpbasis, full_n_radbasfn, n_grid_pt, l, itype, eig, eigv)
      implicit NONE
      class(t_mpbasis), intent(inout)  :: mpbasis
      integer, intent(in)              :: full_n_radbasfn, n_grid_pt, l, itype
      real, intent(in)                 :: eig(:), eigv(:,:)

      integer :: nn, i

      ! reduced number of basis functions
      nn = mpbasis%num_radbasfn(l, itype)

      DO i = 1, n_grid_pt
         mpbasis%radbasfn_mt(i, 1:nn, l, itype) &
            = MATMUL(mpbasis%radbasfn_mt(i, 1:full_n_radbasfn, l, itype), eigv(:,1:nn))/SQRT(eig(:nn))
      END DO
   end subroutine mpbasis_trafo_to_orthonorm_bas

   subroutine mpbasis_add_l0_fun(mpbasis, atoms, hybrid, n_grid_pt, l, itype, gridf)
      use m_types_setup
      use m_types_hybrid
      USE m_intgrf, ONLY: intgrf
      implicit none
      class(t_mpbasis), intent(inout) :: mpbasis
      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                         :: i, j, nn
      REAL, ALLOCATABLE               :: basmhlp(:,:,:,:)
      real                            :: norm

      nn = mpbasis%num_radbasfn(l, itype)
      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

         ! 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
         DO i = nn, 1, -1
            DO j = i + 1, nn + 1
               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)

            END DO

            ! 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
         END DO
         nn = nn + 1
         mpbasis%num_radbasfn(l, itype) = nn
      END IF
   end subroutine mpbasis_add_l0_fun
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442

   subroutine mpbasis_reduce_linear_dep(mpbasis, atoms, mpi, hybrid, l, itype, gridf)
      use m_types_setup
      use m_types_hybrid
      use m_types_mpi
      implicit none
      class(t_mpbasis)              :: mpbasis
      type(t_atoms), intent(in)     :: atoms
      type(t_mpi), intent(in)       :: mpi
      type(t_hybrid), intent(in)    :: hybrid
      integer, intent(in)           :: l, itype

      REAL, ALLOCATABLE             :: olap(:,:), eig(:), eigv(:,:)
      REAL                          :: gridf(:,:)
      integer                       :: full_n_radbasfn, n_grid_pt

      full_n_radbasfn = mpbasis%num_radbasfn(l, itype)
      n_grid_pt = atoms%jri(itype)

      ! Calculate overlap
      call mpbasis%calc_olap_radbasfn(atoms, l, itype, gridf, olap)

      ! Diagonalize
      call mpbasis_diagonialize_olap(olap, eig, eigv)

      call mpbasis%filter_radbasfn(l, itype, full_n_radbasfn, eig, eigv)

      call mpbasis%trafo_to_orthonorm_bas(full_n_radbasfn, n_grid_pt, l, itype, eig, eigv)

      ! Add constant function to l=0 basis and then do a Gram-Schmidt orthonormalization
      call mpbasis%add_l0_fun(atoms, hybrid, n_grid_pt, l, itype, gridf)

      ! Check orthonormality of product basis
      call mpbasis%check_orthonormality(atoms, mpi, l, itype, gridf)

      deallocate(olap, eigv, eig)
   end subroutine
Matthias Redies's avatar
Matthias Redies committed
443
end module m_types_mpbasis