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

      CONTAINS

5 6 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
         SUBROUTINE checkolap(atoms, hybdat,&
        &                  hybrid,&
        &                  nkpti, kpts,&
        &                  dimension, mpi, skip_kpt,&
        &                  input, sym, noco,&
        &                  cell, lapw, jsp)
            USE m_util, ONLY: intgrf, intgrf_init, chr, sphbessel, harmonicsr
            USE m_constants
            USE m_types
            USE m_io_hybrid
            IMPLICIT NONE

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

            TYPE(t_mpi), INTENT(IN)         :: mpi
            TYPE(t_dimension), INTENT(IN)   :: dimension
            TYPE(t_hybrid), INTENT(IN)      :: hybrid
            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

            ! - arrays -
            LOGICAL, INTENT(IN)     ::  skip_kpt(nkpti)

            ! - local scalars -
            INTEGER                 ::  i, itype, iatom, ikpt, ineq, igpt, iband
            INTEGER                 ::  irecl_cmt
            INTEGER                 ::  j, m
            INTEGER                 ::  l
            INTEGER                 :: lm, lm1
            INTEGER                 ::  n, nred, nbasfcn

            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)
            REAL                    ::  bkpt(3)
            REAL                    ::  rarr(maxval(hybrid%nbands))
            REAL                    ::  rtaual(3)
            REAL, ALLOCATABLE   ::  olapcb(:)
            REAL, ALLOCATABLE   :: olapcv_avg(:, :, :, :), olapcv_max(:, :, :, :)
            TYPE(t_mat), ALLOCATABLE :: z(:)

            COMPLEX                 ::  cmt(dimension%neigd, hybrid%maxlmindx, atoms%nat, nkpti)
            COMPLEX                 ::  y((atoms%lmaxd + 1)**2)
            COMPLEX, ALLOCATABLE   ::  olapcv(:, :), olapww(:, :)
            COMPLEX, ALLOCATABLE   ::  carr1(:, :), carr2(:, :), carr3(:, :)

            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'/)
            LOGICAL                 ::  l_mism = .true.

            ALLOCATE (z(nkpti))
            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)
               call z(ikpt)%alloc(sym%invs, nbasfcn, dimension%neigd)
            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)
93 94
            END DO

95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
            IF (mpi%irank == 0) WRITE (6, '(/A)') ' Overlap <core|core>'
            DO itype = 1, atoms%ntype
               IF (atoms%ntype > 1 .AND. mpi%irank == 0) &
            &     WRITE (6, '(A,I3)') ' Atom type', itype
               DO l = 0, hybdat%lmaxc(itype)
                  DO i = 1, hybdat%nindxc(l, itype)
                     IF (mpi%irank == 0)&
              &        WRITE (6, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
                     DO j = 1, i
                        integrand = hybdat%core1(:, i, l, itype)*hybdat%core1(:, j, l, itype)&
               &                  + hybdat%core2(:, i, l, itype)*hybdat%core2(:, j, l, itype)
                        IF (mpi%irank == 0) WRITE (6, '(F10.6)', advance='no')&
               &           intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
                     END DO
                     IF (mpi%irank == 0) WRITE (6, *)
                  END DO
               END DO
112
            END DO
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163

            IF (mpi%irank == 0) WRITE (6, '(/A)') ' Overlap <core|basis>'
            ALLOCATE (olapcb(hybrid%maxindx), olapcv(maxval(hybrid%nbands), nkpti),&
           &          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))

            DO itype = 1, atoms%ntype
               IF (atoms%ntype > 1 .AND. mpi%irank == 0) &
            &     WRITE (6, '(A,I3)') ' Atom type', itype
               DO l = 0, hybdat%lmaxc(itype)
                  IF (l > atoms%lmax(itype)) EXIT ! very improbable case
                  IF (mpi%irank == 0) &
             &        WRITE (6, "(9X,'u(',A,')',4X,'udot(',A,')',:,3X,'ulo(',A,"//&
             &                "') ...')") (lchar(l), i=1, min(3, hybrid%nindx(l, itype)))
                  DO i = 1, hybdat%nindxc(l, itype)
                     IF (mpi%irank == 0)&
              &        WRITE (6, '(1x,I1,A,2X)', advance='no') i + l, lchar(l)
                     DO j = 1, hybrid%nindx(l, itype)

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

                        olapcb(j) = &
               &              intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)

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

                     lm = sum((/(hybrid%nindx(j, itype)*(2*j + 1), j=0, l - 1)/))
                     iatom = sum(atoms%neq(1:itype - 1)) + 1 ! take first of group of equivalent atoms
                     DO m = -l, l
                        olapcv = 0
                        DO j = 1, hybrid%nindx(l, itype)
                           lm = lm + 1
                           olapcv(:, :) = olapcv(:, :) + &
                &                        olapcb(j)*cmt(:maxval(hybrid%nbands), lm, iatom, :nkpti)
                        END DO
                        rdum = sum(abs(olapcv(:, :))**2)
                        rdum1 = maxval(abs(olapcv(:, :)))
                        iarr = maxloc(abs(olapcv(:, :)))
                        olapcv_avg(m, i, l, itype) = &
               &                sqrt(rdum/nkpti/sum(hybrid%nbands(:nkpti))*nkpti)
                        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
164 165
            END DO

166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
            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)') &
               &                                        olapcv_avg(-l:l, i, l, itype)
                     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)//&
               &                 '(F10.6,'' ('',I3.3,''/'',I4.3,'')''))')&
               &                          (olapcv_max(m, i, l, itype),&
               &                           olapcv_loc(:, m, i, l, itype), m=-l, l)
                     END DO
190
                  END DO
191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
               END DO
            END IF ! mpi%irank == 0

            DEALLOCATE (olapcb, olapcv, olapcv_avg, olapcv_max, olapcv_loc)

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

            DO itype = 1, atoms%ntype
               IF (atoms%ntype > 1 .AND. mpi%irank == 0) &
            &     WRITE (6, '(A,I3)') ' Atom type', itype
               DO l = 0, atoms%lmax(itype)
                  DO i = 1, hybrid%nindx(l, itype)
                     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)&
               &                  + hybdat%bas2(:, i, l, itype)*hybdat%bas2(:, j, l, itype)

                        IF (mpi%irank == 0) WRITE (6, '(F10.6)', advance='no')&
               &              intgrf(integrand, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
                     END DO
                     IF (mpi%irank == 0) WRITE (6, *)
                  END DO
               END DO
223
            END DO
224 225 226 227 228 229 230 231 232 233

            IF (.not. l_mism) RETURN

            IF (mpi%irank == 0) WRITE (6, '(/A)') &
           &          'Mismatch of wave functions at the MT-sphere boundaries'
            ALLOCATE (carr1(maxval(hybrid%nbands), (atoms%lmaxd + 1)**2))
            ALLOCATE (carr2(maxval(hybrid%nbands), (atoms%lmaxd + 1)**2))
            ALLOCATE (carr3(maxval(hybrid%nbands), (atoms%lmaxd + 1)**2))
            DO ikpt = 1, nkpti
               call read_z(z(ikpt), ikpt)
234
            END DO
235 236 237 238 239 240 241 242 243 244 245 246 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 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309

            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)
                        gpt(1) = lapw%k1(igpt, jsp)
                        gpt(2) = lapw%k2(igpt, jsp)
                        gpt(3) = lapw%k3(igpt, jsp)

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

                        qnorm = sqrt(sum(q**2))
                        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
                              DO iband = 1, hybrid%nbands(ikpt)
                                 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
                           DO n = 1, hybrid%nindx(l, itype)
                              lm1 = lm1 + 1
                              rdum = hybdat%bas1(atoms%jri(itype), n, l, itype)/atoms%rmt(itype)
                              DO iband = 1, hybrid%nbands(ikpt)
                                 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,'')'')') &
                     !    &              ikpt,sum(rarr(:1)**2/nbands(ikpt)),maxval(rarr(:1))
                     !             CALL writeout(outtext,mpi%irank)
310
!             IF( iatom .eq. 6 ) THEN
311
!               cdum = exp(2*pi*img*dot_product(bkf(:,ikpt),(/0.0,0.0,1.0/) ))
312
!               lm = 0
313 314 315 316 317 318 319 320 321 322 323
!               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),
!      &                                              carr2(iband,lm)/(carr3(iband,lm))
!                   END DO
!                 END DO
!               END DO
!             END IF

324 325 326
                  END DO
               END DO
            END DO
327

328
         END SUBROUTINE checkolap
329 330

      END MODULE m_checkolap