olap.F90 14.2 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 157 158 159 160
   SUBROUTINE wfolap_init(olappw, olapmt, gpt,&
  &                       atoms, hybrid, cell,&
  &                       bas1, bas2)

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

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

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

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

Matthias Redies's avatar
Matthias Redies committed
179
      CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
180
      olapmt = 0
Matthias Redies's avatar
Matthias Redies committed
181 182
      DO itype = 1, atoms%ntype
         DO l = 0, atoms%lmax(itype)
183
            nn = hybrid%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
184 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),&
         &                            atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
191
!               olapmt(n2,n1,l,itype) = olapmt(n1,n2,l,itype)
Matthias Redies's avatar
Matthias Redies committed
192
               END DO
193
            END DO
Matthias Redies's avatar
Matthias Redies committed
194
         END DO
195 196
      END DO

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

   END SUBROUTINE wfolap_init
200

Matthias Redies's avatar
Matthias Redies committed
201
   FUNCTION wfolap_inv(cmt1, cpw1, cmt2, cpw2, ngpt1, ngpt2, olappw, olapmt, atoms, hybrid)
202 203 204 205

      USE m_wrapper
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
206 207
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms
208

Matthias Redies's avatar
Matthias Redies committed
209
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
210
      COMPLEX                :: wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
211 212
      INTEGER, INTENT(IN)     :: ngpt1, ngpt2
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
213 214 215 216
      COMPLEX, INTENT(IN)     :: cmt1(:,:), cmt2(:,:)
      REAL, INTENT(IN)        :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      REAL, INTENT(IN)        :: olappw(:,:)
Matthias Redies's avatar
Matthias Redies committed
217
      REAL, INTENT(IN)        :: olapmt(maxval(hybrid%num_radfun_per_l), maxval(hybrid%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
Daniel Wortmann's avatar
Daniel Wortmann committed
218
!     - local -
Matthias Redies's avatar
Matthias Redies committed
219 220 221
      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
222 223

      wfolap_inv = 0
Matthias Redies's avatar
Matthias Redies committed
224 225 226 227 228 229 230
      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
231
                  nn = hybrid%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
232 233 234 235 236 237
                  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
238
            END DO
Matthias Redies's avatar
Matthias Redies committed
239
         END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
240
      END DO
Matthias Redies's avatar
Matthias Redies committed
241 242

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

244 245
!       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
246
!
247 248
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
Daniel Wortmann's avatar
Daniel Wortmann committed
249 250 251 252
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

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

Daniel Wortmann's avatar
Daniel Wortmann committed
256 257 258
      USE m_wrapper
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
259 260
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
261

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

Daniel Wortmann's avatar
Daniel Wortmann committed
276
      wfolap_noinv = 0
Matthias Redies's avatar
Matthias Redies committed
277 278 279 280 281 282 283
      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
284
                  nn = hybrid%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
285 286 287 288 289 290
                  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
291
            END DO
Matthias Redies's avatar
Matthias Redies committed
292
         END DO
293
      END DO
Matthias Redies's avatar
Matthias Redies committed
294 295

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

297 298
!       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
299
!
300 301
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
302 303 304 305
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
306
   END FUNCTION wfolap_noinv
307

Matthias Redies's avatar
Matthias Redies committed
308 309
   FUNCTION wfolap1(cmt1, cpw1, cmt2, cpw2, ngpt1, ngpt2, olappw, olapmt,&
  &                atoms, hybrid)
310

Matthias Redies's avatar
Matthias Redies committed
311
      USE m_types
312
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
313 314 315 316 317

      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms

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

334
!     - local -
Matthias Redies's avatar
Matthias Redies committed
335 336
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn

337
      wfolap1 = 0
Matthias Redies's avatar
Matthias Redies committed
338 339 340 341 342 343 344
      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
345
                  nn = hybrid%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
346 347 348 349 350 351
                  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
352
            END DO
Matthias Redies's avatar
Matthias Redies committed
353
         END DO
354 355
      END DO

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

Matthias Redies's avatar
Matthias Redies committed
358
   END FUNCTION wfolap1
359

Matthias Redies's avatar
Matthias Redies committed
360 361 362 363
   FUNCTION wfolap2(cmt1, cpw1, cmt2, cpw2, ngpt1, ngpt2, olappw, olapmt,&
  &                atoms, hybrid)
      USE m_types
      IMPLICIT NONE
364

Matthias Redies's avatar
Matthias Redies committed
365 366
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms
367

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

388
      wfolap2 = 0
Matthias Redies's avatar
Matthias Redies committed
389 390 391 392 393 394 395
      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
396
                  nn = hybrid%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
397 398 399 400 401 402
                  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
403
            END DO
Matthias Redies's avatar
Matthias Redies committed
404
         END DO
405 406
      END DO

Matthias Redies's avatar
Matthias Redies committed
407 408 409 410
      wfolap2 = wfolap2 + dot_product(cpw1, matmul(olappw, cpw2))

   END FUNCTION wfolap2
END MODULE m_olap