olap.F90 14.4 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1
MODULE m_olap
2

Matthias Redies's avatar
Matthias Redies committed
3
CONTAINS
4 5 6 7 8

!     Calculates plane-wave overlap matrix olap defined by GPT(1:3,1:NGPT).
!     (Muffin-tin spheres are cut out.)
!     olap_pw calculates full overlap matrix

Matthias Redies's avatar
Matthias Redies committed
9 10 11
   SUBROUTINE olap_pw(olap, gpt, ngpt, atoms, cell)

      USE m_constants
12 13
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
14 15
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms
16 17

!     - scalars -
Matthias Redies's avatar
Matthias Redies committed
18
      INTEGER, INTENT(IN)       :: ngpt
19
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
20
      INTEGER, INTENT(IN)       :: gpt(:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
21
      TYPE(t_mat)              :: olap
22
!     - local -
Matthias Redies's avatar
Matthias Redies committed
23 24 25
      INTEGER                  :: i, j, itype, icent, ineq
      REAL                     :: g, r, fgr
      COMPLEX, PARAMETER        :: img = (0.0, 1.0)
26
      INTEGER                  :: dg(3)
Matthias Redies's avatar
Matthias Redies committed
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

      DO i = 1, ngpt
         DO j = 1, i
            dg = gpt(:, j) - gpt(:, i)
            g = gptnorm(dg, cell%bmat)
            IF (g == 0) THEN
               DO itype = 1, atoms%ntype
                  r = atoms%rmt(itype)
                  if (olap%l_real) THEN
                     olap%data_r(i, j) = olap%data_r(i, j) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
                  else
                     olap%data_c(i, j) = olap%data_c(i, j) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
                  endif
               END DO
            ELSE
               icent = 0
               DO itype = 1, atoms%ntype
                  r = g*atoms%rmt(itype)
                  fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
                  DO ineq = 1, atoms%neq(itype)
                     icent = icent + 1
                     if (olap%l_real) THEN
                        olap%data_r(i, j) = olap%data_r(i, j) - fgr*exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
                     else
                        olap%data_c(i, j) = olap%data_c(i, j) - fgr*exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
                     endif
                  END DO
               END DO
            END IF
            if (olap%l_real) THEN
               IF (i == j) olap%data_r(i, j) = olap%data_r(i, j) + 1
               olap%data_r(j, i) = olap%data_r(i, j)
            else
               IF (i == j) olap%data_c(i, j) = olap%data_c(i, j) + 1
               olap%data_c(j, i) = conjg(olap%data_c(i, j))
            endif
         END DO
64 65
      END DO

Matthias Redies's avatar
Matthias Redies committed
66
   END SUBROUTINE olap_pw
67

Matthias Redies's avatar
Matthias Redies committed
68
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69 70 71

!     olap_pwp  calculates upper triangular part of overlap matrix

Matthias Redies's avatar
Matthias Redies committed
72
   SUBROUTINE olap_pwp(l_real, olap_r, olap_c, gpt, ngpt, atoms, cell)
73

Matthias Redies's avatar
Matthias Redies committed
74
      USE m_constants
75 76
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
77 78
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms
79 80

!     - scalars -
Matthias Redies's avatar
Matthias Redies committed
81
      INTEGER, INTENT(IN)       :: ngpt
82
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
83
      INTEGER, INTENT(IN)       :: gpt(:,:)
84

Matthias Redies's avatar
Matthias Redies committed
85 86 87
      LOGICAL, INTENT(IN)       :: l_real
      REAL, INTENT(OUT)         ::  olap_r(ngpt*(ngpt + 1)/2)
      COMPLEX, INTENT(OUT)      ::  olap_c(ngpt*(ngpt + 1)/2)
88
!     - local -
Matthias Redies's avatar
Matthias Redies committed
89 90 91
      INTEGER                  :: i, j, k, itype, icent, ineq
      REAL                     :: g, r, fgr
      COMPLEX, PARAMETER        :: img = (0.0, 1.0)
92 93
      INTEGER                  :: dg(3)

Daniel Wortmann's avatar
Daniel Wortmann committed
94
      if (l_real) THEN
Matthias Redies's avatar
Matthias Redies committed
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
         k = 0
         DO i = 1, ngpt
            DO j = 1, i
               k = k + 1
               dg = gpt(:, i) - gpt(:, j)
               g = gptnorm(dg, cell%bmat)
               olap_r(k) = 0
               IF (g == 0) THEN
                  DO itype = 1, atoms%ntype
                     r = atoms%rmt(itype)
                     olap_r(k) = olap_r(k) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
                  END DO
               ELSE
                  icent = 0
                  DO itype = 1, atoms%ntype
                     r = g*atoms%rmt(itype)
                     fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
                     DO ineq = 1, atoms%neq(itype)
                        icent = icent + 1
                        olap_r(k) = olap_r(k) - fgr* &
             &               exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
                     END DO
                  END DO
               END IF
               IF (i == j) olap_r(k) = olap_r(k) + 1
120
            END DO
Matthias Redies's avatar
Matthias Redies committed
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
         END DO
      else
         k = 0
         DO i = 1, ngpt
            DO j = 1, i
               k = k + 1
               dg = gpt(:, i) - gpt(:, j)
               g = gptnorm(dg, cell%bmat)
               olap_c(k) = 0
               IF (g == 0) THEN
                  DO itype = 1, atoms%ntype
                     r = atoms%rmt(itype)
                     olap_c(k) = olap_c(k) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
                  END DO
               ELSE
                  icent = 0
                  DO itype = 1, atoms%ntype
                     r = g*atoms%rmt(itype)
                     fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
                     DO ineq = 1, atoms%neq(itype)
                        icent = icent + 1
                        olap_c(k) = olap_c(k) - fgr* &
             &               exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
                     END DO
                  END DO
               END IF
               IF (i == j) olap_c(k) = olap_c(k) + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
148
            END DO
Matthias Redies's avatar
Matthias Redies committed
149
         END DO
150

Matthias Redies's avatar
Matthias Redies committed
151 152
      endif
   END SUBROUTINE olap_pwp
153 154 155

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Matthias Redies's avatar
Matthias Redies committed
156
   SUBROUTINE wfolap_init(olappw, olapmt, gpt,&
157
  &                       atoms, mpbasis, hybrid, cell,&
Matthias Redies's avatar
Matthias Redies committed
158 159
  &                       bas1, bas2)

Matthias Redies's avatar
Matthias Redies committed
160
      USE m_intgrf, ONLY: intgrf, intgrf_init
161 162
      USE m_types
      IMPLICIT NONE
163
      TYPE(t_mpbasis), intent(in) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
164 165 166 167
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms

168
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
169
      INTEGER, INTENT(IN)       :: gpt(:, :)!(3,ngpt)
170 171 172
      REAL, INTENT(IN)         ::  bas1(atoms%jmtd, maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype),&
     &                            bas2(atoms%jmtd, maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
      REAL, INTENT(OUT)         :: olapmt(maxval(mpbasis%num_radfun_per_l), maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
Matthias Redies's avatar
Matthias Redies committed
173
      TYPE(t_mat), INTENT(INOUT):: olappw
174 175

!     - local -
Matthias Redies's avatar
Matthias Redies committed
176
      INTEGER                  :: itype, l, nn, n1, n2
177

Matthias Redies's avatar
Matthias Redies committed
178
      REAL, ALLOCATABLE         :: gridf(:, :)
179

Matthias Redies's avatar
Matthias Redies committed
180
      CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
181
      olapmt = 0
Matthias Redies's avatar
Matthias Redies committed
182 183
      DO itype = 1, atoms%ntype
         DO l = 0, atoms%lmax(itype)
184
            nn = mpbasis%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
185 186 187 188 189 190
            DO n2 = 1, nn
               DO n1 = 1, nn!n2
                  !IF( n1 .gt. 2 .or. n2 .gt. 2) CYCLE
                  olapmt(n1, n2, l, itype) = intgrf( &
         &                            bas1(:, n1, l, itype)*bas1(:, n2, l, itype)&
         &                           + bas2(:, n1, l, itype)*bas2(:, n2, l, itype),&
Matthias Redies's avatar
Matthias Redies committed
191
         &                            atoms, itype, gridf)
192
!               olapmt(n2,n1,l,itype) = olapmt(n1,n2,l,itype)
Matthias Redies's avatar
Matthias Redies committed
193
               END DO
194
            END DO
Matthias Redies's avatar
Matthias Redies committed
195
         END DO
196 197
      END DO

Matthias Redies's avatar
Matthias Redies committed
198 199 200
      CALL olap_pw(olappw, gpt, size(gpt, 2), atoms, cell)

   END SUBROUTINE wfolap_init
201

202
   FUNCTION wfolap_inv(cmt1, cpw1, cmt2, cpw2, ngpt1, ngpt2, olappw, olapmt, atoms, mpbasis,hybrid)
203 204 205 206

      USE m_wrapper
      USE m_types
      IMPLICIT NONE
207
      TYPE(t_mpbasis), intent(in) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
208 209
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms
210

Matthias Redies's avatar
Matthias Redies committed
211
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
212
      COMPLEX                :: wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
213 214
      INTEGER, INTENT(IN)     :: ngpt1, ngpt2
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
215 216 217 218
      COMPLEX, INTENT(IN)     :: cmt1(:,:), cmt2(:,:)
      REAL, INTENT(IN)        :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      REAL, INTENT(IN)        :: olappw(:,:)
219
      REAL, INTENT(IN)        :: olapmt(maxval(mpbasis%num_radfun_per_l), maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
Daniel Wortmann's avatar
Daniel Wortmann committed
220
!     - local -
Matthias Redies's avatar
Matthias Redies committed
221 222 223
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
      COMPLEX                :: carr(ngpt1), cdum
      REAL                   :: rarr1(ngpt1), rarr2(ngpt1), rdum1, rdum2
Daniel Wortmann's avatar
Daniel Wortmann committed
224 225

      wfolap_inv = 0
Matthias Redies's avatar
Matthias Redies committed
226 227 228 229 230 231 232
      iatom = 0
      DO itype = 1, atoms%ntype
         DO ieq = 1, atoms%neq(itype)
            iatom = iatom + 1
            lm = 0
            DO l = 0, atoms%lmax(itype)
               DO M = -l, l
233
                  nn = mpbasis%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
234 235 236 237 238 239
                  wfolap_inv = wfolap_inv + &
         &                 dot_product(cmt1(lm + 1:lm + nn, iatom),&
         &                               matmul(olapmt(:nn, :nn, l, itype),&
         &                                       cmt2(lm + 1:lm + nn, iatom)))
                  lm = lm + nn
               END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
240
            END DO
Matthias Redies's avatar
Matthias Redies committed
241
         END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
242
      END DO
Matthias Redies's avatar
Matthias Redies committed
243 244

      wfolap_inv = wfolap_inv + dot_product(cpw1, matmul(olappw, cpw2))
Daniel Wortmann's avatar
Daniel Wortmann committed
245

246 247
!       CALL dgemv('N',ngpt1,ngpt2,1.0,olappw,ngpt1,real(cpw2),1,0.0,rarr1,1)
!       CALL dgemv('N',ngpt1,ngpt2,1.0,olappw,ngpt1,aimag(cpw2),1,0.0,rarr2,1)
Matthias Redies's avatar
Matthias Redies committed
248
!
249 250
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
Daniel Wortmann's avatar
Daniel Wortmann committed
251 252 253 254
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
255
   END FUNCTION wfolap_inv
256
   FUNCTION wfolap_noinv(cmt1, cpw1, cmt2, cpw2, ngpt1, ngpt2, olappw, olapmt, atoms, mpbasis, hybrid)
Matthias Redies's avatar
Matthias Redies committed
257

Daniel Wortmann's avatar
Daniel Wortmann committed
258 259 260
      USE m_wrapper
      USE m_types
      IMPLICIT NONE
261
      TYPE(t_mpbasis), intent(in) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
262 263
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
264

Matthias Redies's avatar
Matthias Redies committed
265
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
266
      COMPLEX                :: wfolap_noinv
Matthias Redies's avatar
Matthias Redies committed
267 268
      INTEGER, INTENT(IN)     :: ngpt1, ngpt2
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
269 270 271 272
      COMPLEX, INTENT(IN)     :: cmt1(:,:),cmt2(:,:)
      COMPLEX, INTENT(IN)     :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      COMPLEX, INTENT(IN)     :: olappw(:,:)
273
      REAL, INTENT(IN)        :: olapmt(maxval(mpbasis%num_radfun_per_l), maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
274
!     - local -
Matthias Redies's avatar
Matthias Redies committed
275 276 277
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
      COMPLEX                :: carr(ngpt1), cdum
      REAL                   :: rarr1(ngpt1), rarr2(ngpt1), rdum1, rdum2
278

Daniel Wortmann's avatar
Daniel Wortmann committed
279
      wfolap_noinv = 0
Matthias Redies's avatar
Matthias Redies committed
280 281 282 283 284 285 286
      iatom = 0
      DO itype = 1, atoms%ntype
         DO ieq = 1, atoms%neq(itype)
            iatom = iatom + 1
            lm = 0
            DO l = 0, atoms%lmax(itype)
               DO M = -l, l
287
                  nn = mpbasis%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
288 289 290 291 292 293
                  wfolap_noinv = wfolap_noinv + &
         &                 dot_product(cmt1(lm + 1:lm + nn, iatom),&
         &                               matmul(olapmt(:nn, :nn, l, itype),&
         &                                       cmt2(lm + 1:lm + nn, iatom)))
                  lm = lm + nn
               END DO
294
            END DO
Matthias Redies's avatar
Matthias Redies committed
295
         END DO
296
      END DO
Matthias Redies's avatar
Matthias Redies committed
297 298

      wfolap_noinv = wfolap_noinv + dot_product(cpw1, matmul(olappw, cpw2))
299

300 301
!       CALL dgemv('N',ngpt1,ngpt2,1.0,olappw,ngpt1,real(cpw2),1,0.0,rarr1,1)
!       CALL dgemv('N',ngpt1,ngpt2,1.0,olappw,ngpt1,aimag(cpw2),1,0.0,rarr2,1)
Matthias Redies's avatar
Matthias Redies committed
302
!
303 304
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
305 306 307 308
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
309
   END FUNCTION wfolap_noinv
310

Matthias Redies's avatar
Matthias Redies committed
311
   FUNCTION wfolap1(cmt1, cpw1, cmt2, cpw2, ngpt1, ngpt2, olappw, olapmt,&
312
  &                atoms, mpbasis,hybrid)
313

Matthias Redies's avatar
Matthias Redies committed
314
      USE m_types
315
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
316

317
      TYPE(t_mpbasis), intent(in) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
318 319 320 321
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms

!     -scalars -
322
      COMPLEX                :: wfolap1
Matthias Redies's avatar
Matthias Redies committed
323
      INTEGER, INTENT(IN)     :: ngpt1, ngpt2
324
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
325
      COMPLEX, INTENT(IN)     :: cmt1(:,:), cmt2(:,:)
326
#if ( defined(CPP_INVERSION) )
Matthias Redies's avatar
Matthias Redies committed
327
      REAL, INTENT(IN)        :: cpw1(:), cpw2(:)
328
#else
Matthias Redies's avatar
Matthias Redies committed
329
      COMPLEX, INTENT(IN)     :: cpw1(:), cpw2(:)
330 331
#endif
#if ( defined(CPP_INVERSION) )
Matthias Redies's avatar
Matthias Redies committed
332
      REAL, INTENT(IN)        :: olappw(:,:)
333
#else
Matthias Redies's avatar
Matthias Redies committed
334
      COMPLEX, INTENT(IN)     :: olappw(:,:)
335
#endif
336
      REAL, INTENT(IN)        :: olapmt(maxval(mpbasis%num_radfun_per_l), maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
Matthias Redies's avatar
Matthias Redies committed
337

338
!     - local -
Matthias Redies's avatar
Matthias Redies committed
339 340
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn

341
      wfolap1 = 0
Matthias Redies's avatar
Matthias Redies committed
342 343 344 345 346 347 348
      iatom = 0
      DO itype = 1, atoms%ntype
         DO ieq = 1, atoms%neq(itype)
            iatom = iatom + 1
            lm = 0
            DO l = 0, atoms%lmax(itype)
               DO M = -l, l
349
                  nn = mpbasis%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
350 351 352 353 354 355
                  wfolap1 = wfolap1 + &
         &                  dot_product(cmt1(lm + 1:lm + nn, iatom),&
         &                                matmul(olapmt(:nn, :nn, l, itype),&
         &                                        cmt2(lm + 1:lm + nn, iatom)))
                  lm = lm + nn
               END DO
356
            END DO
Matthias Redies's avatar
Matthias Redies committed
357
         END DO
358 359
      END DO

Matthias Redies's avatar
Matthias Redies committed
360
      wfolap1 = wfolap1 + dot_product(cpw1, matmul(olappw, cpw2))
361

Matthias Redies's avatar
Matthias Redies committed
362
   END FUNCTION wfolap1
363

Matthias Redies's avatar
Matthias Redies committed
364
   FUNCTION wfolap2(cmt1, cpw1, cmt2, cpw2, ngpt1, ngpt2, olappw, olapmt,&
365
  &                atoms, mpbasis, hybrid)
Matthias Redies's avatar
Matthias Redies committed
366 367
      USE m_types
      IMPLICIT NONE
368

369
      TYPE(t_mpbasis), intent(in) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
370 371
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms
372

Matthias Redies's avatar
Matthias Redies committed
373
!     - scalars -
374
      COMPLEX                :: wfolap2
Matthias Redies's avatar
Matthias Redies committed
375 376
      INTEGER, INTENT(IN)     :: ngpt1, ngpt2
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
377
      COMPLEX, INTENT(IN)     :: cmt1(:,:),cmt2(:,:)
378
! #if ( defined(CPP_INVERSION) )
Matthias Redies's avatar
Matthias Redies committed
379
!       REAL,INTENT(IN)        :: cpw1(:)
380
! #else
Matthias Redies's avatar
Matthias Redies committed
381
      COMPLEX, INTENT(IN)     :: cpw1(:)
382
! #endif
Matthias Redies's avatar
Matthias Redies committed
383
      COMPLEX, INTENT(IN)     :: cpw2(:)
384
#if ( defined(CPP_INVERSION) )
Matthias Redies's avatar
Matthias Redies committed
385
      REAL, INTENT(IN)        :: olappw(:,:)
386
#else
Matthias Redies's avatar
Matthias Redies committed
387
      COMPLEX, INTENT(IN)     :: olappw(:,:)
388
#endif
389
      REAL, INTENT(IN)        :: olapmt(maxval(mpbasis%num_radfun_per_l), maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
390
!     - local -
Matthias Redies's avatar
Matthias Redies committed
391 392
      INTEGER                :: itype, ieq, ic, l, m, lm, nn

393
      wfolap2 = 0
Matthias Redies's avatar
Matthias Redies committed
394 395 396 397 398 399 400
      ic = 0
      DO itype = 1, atoms%ntype
         DO ieq = 1, atoms%neq(itype)
            ic = ic + 1
            lm = 0
            DO l = 0, atoms%lmax(itype)
               DO M = -l, l
401
                  nn = mpbasis%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
402 403 404 405 406 407
                  wfolap2 = wfolap2 + &
         &                 dot_product(cmt1(lm + 1:lm + nn, ic),&
         &                               matmul(olapmt(:nn, :nn, l, itype),&
         &                                       cmt2(lm + 1:lm + nn, ic)))
                  lm = lm + nn
               END DO
408
            END DO
Matthias Redies's avatar
Matthias Redies committed
409
         END DO
410 411
      END DO

Matthias Redies's avatar
Matthias Redies committed
412 413 414 415
      wfolap2 = wfolap2 + dot_product(cpw1, matmul(olappw, cpw2))

   END FUNCTION wfolap2
END MODULE m_olap