checkolap.F90 14.6 KB
Newer Older
1 2 3 4
      MODULE m_checkolap

      CONTAINS

5
         SUBROUTINE checkolap(atoms, hybdat,&
6
                           mpdata,hybinp,&
7
                           nkpti, kpts,&
8
                            mpi, &
9 10
                           input, sym, noco,&
                           cell, lapw, jsp)
11 12
            USE m_util, ONLY: chr, sphbessel, harmonicsr
            use m_intgrf, only:  intgrf, intgrf_init
13 14
            USE m_constants
            USE m_types
15
            USE m_io_hybinp
16 17
            USE m_types_hybdat

18 19 20 21 22
            IMPLICIT NONE

            TYPE(t_hybdat), INTENT(IN)   :: hybdat

            TYPE(t_mpi), INTENT(IN)         :: mpi
23
            TYPE(t_mpdata), intent(in) :: mpdata
24
            TYPE(t_hybinp), INTENT(IN)      :: hybinp
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
            TYPE(t_input), INTENT(IN)       :: input
            TYPE(t_noco), INTENT(IN)        :: noco
            TYPE(t_sym), INTENT(IN)         :: sym
            TYPE(t_cell), INTENT(IN)        :: cell
            TYPE(t_kpts), INTENT(IN)        :: kpts
            TYPE(t_atoms), INTENT(IN)       :: atoms
            TYPE(t_lapw), INTENT(INOUT)     :: lapw

            ! - scalars -
            INTEGER, INTENT(IN)     :: jsp
            INTEGER, INTENT(IN)     ::  nkpti


            ! - local scalars -
            INTEGER                 ::  i, itype, iatom, ikpt, ineq, igpt, iband
            INTEGER                 ::  j, m
            INTEGER                 ::  l
            INTEGER                 :: lm, lm1
Matthias Redies's avatar
Matthias Redies committed
43
            INTEGER                 ::  n, nbasfcn
44 45 46 47 48 49 50 51 52 53 54 55 56 57

            REAL                    ::  rdum, rdum1
            REAL                    ::  qnorm

            COMPLEX                 ::  cexp, cdum
            COMPLEX, PARAMETER     ::  img = (0.0, 1.0)

            ! -local arrays -
            INTEGER                 ::  iarr(2), gpt(3)
            INTEGER, ALLOCATABLE   ::  olapcv_loc(:, :, :, :, :)

            REAL                    ::  sphbes(0:atoms%lmaxd)
            REAL                    ::  q(3)
            REAL                    ::  integrand(atoms%jmtd)
58
            REAL                    ::  rarr(maxval(hybdat%nbands))
59 60 61 62
            REAL, ALLOCATABLE   ::  olapcb(:)
            REAL, ALLOCATABLE   :: olapcv_avg(:, :, :, :), olapcv_max(:, :, :, :)
            TYPE(t_mat), ALLOCATABLE :: z(:)

63
            COMPLEX                 ::  cmt(input%neig, hybdat%maxlmindx, atoms%nat, nkpti)
64
            COMPLEX                 ::  y((atoms%lmaxd + 1)**2)
Matthias Redies's avatar
Matthias Redies committed
65
            COMPLEX, ALLOCATABLE   ::  olapcv(:, :)
66 67 68
            COMPLEX, ALLOCATABLE   ::  carr1(:, :), carr2(:, :), carr3(:, :)

            CHARACTER, PARAMETER    ::  lchar(0:38) =&
69 70 71
                      (/'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'/)
72 73
            LOGICAL                 ::  l_mism = .true.

74
            allocate(z(nkpti))
75 76 77
            DO ikpt = 1, nkpti
               CALL lapw%init(input, noco, kpts, atoms, sym, ikpt, cell, sym%zrfs)
               nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
78
               call z(ikpt)%alloc(sym%invs, nbasfcn, input%neig)
79 80 81 82 83 84 85 86 87 88 89 90
            ENDDO

            IF (mpi%irank == 0) WRITE (6, '(//A)') '### checkolap ###'

            cmt = 0

            ! initialize gridf -> was done in eigen_HF_init
            !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)

            ! read in cmt
            DO ikpt = 1, nkpti
               call read_cmt(cmt(:, :, :, ikpt), ikpt)
91 92
            END DO

93 94 95
            IF (mpi%irank == 0) WRITE (6, '(/A)') ' Overlap <core|core>'
            DO itype = 1, atoms%ntype
               IF (atoms%ntype > 1 .AND. mpi%irank == 0) &
96
                  WRITE (6, '(A,I3)') ' Atom type', itype
97 98 99
               DO l = 0, hybdat%lmaxc(itype)
                  DO i = 1, hybdat%nindxc(l, itype)
                     IF (mpi%irank == 0)&
100
                       WRITE (6, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
101 102
                     DO j = 1, i
                        integrand = hybdat%core1(:, i, l, itype)*hybdat%core1(:, j, l, itype)&
103
                                  + hybdat%core2(:, i, l, itype)*hybdat%core2(:, j, l, itype)
104
                        IF (mpi%irank == 0) WRITE (6, '(F10.6)', advance='no')&
105
                           intgrf(integrand, atoms, itype, hybdat%gridf)
106 107 108 109
                     END DO
                     IF (mpi%irank == 0) WRITE (6, *)
                  END DO
               END DO
110
            END DO
111 112

            IF (mpi%irank == 0) WRITE (6, '(/A)') ' Overlap <core|basis>'
113
            allocate(olapcb(maxval(mpdata%num_radfun_per_l)), olapcv(maxval(hybdat%nbands), nkpti),&
114 115 116
                      olapcv_avg(-hybdat%lmaxcd:hybdat%lmaxcd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype),&
                      olapcv_max(-hybdat%lmaxcd:hybdat%lmaxcd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype),&
                      olapcv_loc(2, -hybdat%lmaxcd:hybdat%lmaxcd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype))
117 118 119

            DO itype = 1, atoms%ntype
               IF (atoms%ntype > 1 .AND. mpi%irank == 0) &
120
                  WRITE (6, '(A,I3)') ' Atom type', itype
121 122 123
               DO l = 0, hybdat%lmaxc(itype)
                  IF (l > atoms%lmax(itype)) EXIT ! very improbable case
                  IF (mpi%irank == 0) &
124
                      WRITE (6, "(9X,'u(',A,')',4X,'udot(',A,')',:,3X,'ulo(',A,"//&
125
                              "') ...')") (lchar(l), i=1, min(3, mpdata%num_radfun_per_l(l, itype)))
126 127
                  DO i = 1, hybdat%nindxc(l, itype)
                     IF (mpi%irank == 0)&
128
                       WRITE (6, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
129
                     DO j = 1, mpdata%num_radfun_per_l(l, itype)
130 131

                        integrand = hybdat%core1(:, i, l, itype)*hybdat%bas1(:, j, l, itype)&
132
                                  + hybdat%core2(:, i, l, itype)*hybdat%bas2(:, j, l, itype)
133 134

                        olapcb(j) = &
135
                              intgrf(integrand, atoms, itype, hybdat%gridf)
136 137

                        IF (mpi%irank == 0)&
138
                          WRITE (6, '(F10.6)', advance='no') olapcb(j)
139 140
                     END DO

141
                     lm = sum([(mpdata%num_radfun_per_l(j, itype)*(2*j + 1), j=0, l - 1)])
142 143 144
                     iatom = sum(atoms%neq(1:itype - 1)) + 1 ! take first of group of equivalent atoms
                     DO m = -l, l
                        olapcv = 0
145
                        DO j = 1, mpdata%num_radfun_per_l(l, itype)
146 147
                           lm = lm + 1
                           olapcv(:, :) = olapcv(:, :) + &
148
                                         olapcb(j)*cmt(:maxval(hybdat%nbands), lm, iatom, :nkpti)
149 150 151 152 153
                        END DO
                        rdum = sum(abs(olapcv(:, :))**2)
                        rdum1 = maxval(abs(olapcv(:, :)))
                        iarr = maxloc(abs(olapcv(:, :)))
                        olapcv_avg(m, i, l, itype) = &
154
                                sqrt(rdum/nkpti/sum(hybdat%nbands(:nkpti))*nkpti)
155 156 157 158 159 160 161
                        olapcv_max(m, i, l, itype) = rdum1
                        olapcv_loc(:, m, i, l, itype) = iarr
                     END DO
                     IF (mpi%irank == 0) WRITE (6, *)

                  END DO
               END DO
162 163
            END DO

164 165 166 167 168 169 170 171
            IF (mpi%irank == 0) THEN
               WRITE (6, '(/A)') ' Average overlap <core|val>'
               DO itype = 1, atoms%ntype
                  IF (atoms%ntype > 1) write (6, '(A,I3)') ' Atom type', itype
                  DO l = 0, hybdat%lmaxc(itype)
                     DO i = 1, hybdat%nindxc(l, itype)
                        WRITE (6, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
                        WRITE (6, '('//chr(2*l + 1)//'F10.6)') &
172
                                                        olapcv_avg(-l:l, i, l, itype)
173 174 175 176 177 178 179 180 181 182 183
                     END DO
                  END DO
               END DO

               WRITE (6, '(/A)') ' Maximum overlap <core|val> at (band/kpoint)'
               DO itype = 1, atoms%ntype
                  IF (atoms%ntype > 1) write (6, '(A,I3)') ' Atom type', itype
                  DO l = 0, hybdat%lmaxc(itype)
                     DO i = 1, hybdat%nindxc(l, itype)
                        WRITE (6, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
                        WRITE (6, '('//chr(2*l + 1)//&
184 185 186
                                 '(F10.6,'' ('',I3.3,''/'',I4.3,'')''))')&
                                          (olapcv_max(m, i, l, itype),&
                                           olapcv_loc(:, m, i, l, itype), m=-l, l)
187
                     END DO
188
                  END DO
189 190 191
               END DO
            END IF ! mpi%irank == 0

192
            deallocate(olapcb, olapcv, olapcv_avg, olapcv_max, olapcv_loc)
193 194 195 196 197

            IF (mpi%irank == 0) WRITE (6, '(/A)') ' Overlap <basis|basis>'

            DO itype = 1, atoms%ntype
               IF (atoms%ntype > 1 .AND. mpi%irank == 0) &
198
                  WRITE (6, '(A,I3)') ' Atom type', itype
199
               DO l = 0, atoms%lmax(itype)
200
                  DO i = 1, mpdata%num_radfun_per_l(l, itype)
201 202 203 204 205 206 207 208 209 210 211 212
                     IF (mpi%irank == 0) THEN
                        SELECT CASE (i)
                        CASE (1)
                           WRITE (6, '(1x,''   u('',A,'')'')', advance='no') lchar(l)
                        CASE (2)
                           WRITE (6, '(1x,''udot('',A,'')'')', advance='no') lchar(l)
                        CASE DEFAULT
                           WRITE (6, '(1x,'' ulo('',A,'')'')', advance='no') lchar(l)
                        END SELECT
                     END IF
                     DO j = 1, i
                        integrand = hybdat%bas1(:, i, l, itype)*hybdat%bas1(:, j, l, itype)&
213
                                  + hybdat%bas2(:, i, l, itype)*hybdat%bas2(:, j, l, itype)
214 215

                        IF (mpi%irank == 0) WRITE (6, '(F10.6)', advance='no')&
216
                              intgrf(integrand, atoms, itype, hybdat%gridf)
217 218 219 220
                     END DO
                     IF (mpi%irank == 0) WRITE (6, *)
                  END DO
               END DO
221
            END DO
222 223 224 225

            IF (.not. l_mism) RETURN

            IF (mpi%irank == 0) WRITE (6, '(/A)') &
226
                      'Mismatch of wave functions at the MT-sphere boundaries'
227 228 229
            allocate(carr1(maxval(hybdat%nbands), (atoms%lmaxd + 1)**2))
            allocate(carr2(maxval(hybdat%nbands), (atoms%lmaxd + 1)**2))
            allocate(carr3(maxval(hybdat%nbands), (atoms%lmaxd + 1)**2))
230
            DO ikpt = 1, nkpti
231
               call read_z(z(ikpt), kpts%nkptf*(jsp - 1) + ikpt)
232
            END DO
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250

            iatom = 0
            DO itype = 1, atoms%ntype
               DO ineq = 1, atoms%neq(itype)
                  iatom = iatom + 1
                  IF (mpi%irank == 0) THEN
                     if (atoms%nat > 1) WRITE (6, '(2X,A,I3)') 'Atom', iatom
                     WRITE (6, '(2X,A)') 'k-point    average      (   maximum    )'
                  END IF

                  DO ikpt = 1, nkpti
                     carr1 = 0; carr2 = 0; carr3 = 0

                     ! calculate k1,k2,k3
                     CALL lapw%init(input, noco, kpts, atoms, sym, ikpt, cell, sym%zrfs)

                     ! PW part
                     DO igpt = 1, lapw%nv(jsp)
Matthias Redies's avatar
Matthias Redies committed
251
                        gpt = lapw%gvec(:,igpt, jsp)
252 253

                        cexp = exp(img*2*pi_const* &
254
                                   dot_product(kpts%bkf(:, ikpt) + gpt, atoms%taual(:, iatom)))
255 256
                        q = matmul(kpts%bkf(:, ikpt) + gpt, cell%bmat)

257
                        qnorm = norm2(q)
258 259 260 261 262 263 264 265
                        call sphbessel(sphbes, atoms%rmt(itype)*qnorm, atoms%lmax(itype))
                        call harmonicsr(y, q, atoms%lmax(itype))
                        y = conjg(y)
                        lm = 0
                        DO l = 0, atoms%lmax(itype)
                           cdum = 4*pi_const*img**l/sqrt(cell%omtil)*sphbes(l)*cexp
                           DO m = -l, l
                              lm = lm + 1
266
                              DO iband = 1, hybdat%nbands(ikpt)
267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282
                                 if (z(1)%l_real) THEN
                                    carr2(iband, lm) = carr2(iband, lm) + cdum*z(ikpt)%data_r(igpt, iband)*y(lm)
                                 Else
                                    carr2(iband, lm) = carr2(iband, lm) + cdum*z(ikpt)%data_c(igpt, iband)*y(lm)
                                 END if
                              end DO
                           END DO
                        END DO
                     END DO

                     ! MT
                     lm = 0
                     lm1 = 0
                     DO l = 0, atoms%lmax(itype)
                        DO m = -l, l
                           lm = lm + 1
283
                           DO n = 1, mpdata%num_radfun_per_l(l, itype)
284 285
                              lm1 = lm1 + 1
                              rdum = hybdat%bas1(atoms%jri(itype), n, l, itype)/atoms%rmt(itype)
286
                              DO iband = 1, hybdat%nbands(ikpt)
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303
                                 carr3(iband, lm) = carr3(iband, lm) + cmt(iband, lm1, iatom, ikpt)*rdum
                              END DO
                           END DO
                        END DO
                     END DO
                     carr1 = carr2 - carr3

                     rarr = 0
                     lm = 0
                     DO l = 0, atoms%lmax(itype)
                        DO m = -l, l
                           lm = lm + 1
                           rarr = rarr + abs(carr1(:, lm))**2
                        END DO
                     END DO
                     rarr = sqrt(rarr/(4*pi_const))
                     !             WRITE(outtext,'(I6,4X,F14.12,''  ('',F14.12,'')'')') &
304
                     !                   ikpt,sum(rarr(:1)**2/nbands(ikpt)),maxval(rarr(:1))
305
                     !             CALL writeout(outtext,mpi%irank)
306
!             IF( iatom .eq. 6 ) THEN
Matthias Redies's avatar
Matthias Redies committed
307
!               cdum = exp(2*pi*img*dot_product(bkf(:,ikpt),[0.0,0.0,1.0] ))
308
!               lm = 0
309 310 311 312 313
!               DO l = 0,lmax(itype)
!                 DO m = -l,l
!                   lm = lm + 1
!                   DO iband = 1,nbands(ikpt)
!                     WRITE(700+ikpt,'(3i4,6f15.10)') iband,l,m,carr2(iband,lm),carr3(iband,lm),
314
!                                                     carr2(iband,lm)/(carr3(iband,lm))
315 316 317 318 319
!                   END DO
!                 END DO
!               END DO
!             END IF

320 321 322
                  END DO
               END DO
            END DO
323

324
         END SUBROUTINE checkolap
325 326

      END MODULE m_checkolap