olap.F90 14.3 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(3, ngpt)
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(3, ngpt)
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
169
170
171
172
      INTEGER, INTENT(IN)       :: gpt(:, :)!(3,ngpt)
      REAL, INTENT(IN)         ::  bas1(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype),&
     &                            bas2(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype)
      REAL, INTENT(OUT)         :: olapmt(hybrid%maxindx, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype)
      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
183
184
185
186
187
188
189
190
      DO itype = 1, atoms%ntype
         DO l = 0, atoms%lmax(itype)
            nn = hybrid%nindx(l, itype)
            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
213
214
215
216
217
218
      INTEGER, INTENT(IN)     :: ngpt1, ngpt2
!     - arrays -
      COMPLEX, INTENT(IN)     :: cmt1(hybrid%maxlmindx, atoms%nat),&
     &                          cmt2(hybrid%maxlmindx, atoms%nat)
      REAL, INTENT(IN)        :: cpw1(ngpt1)
      COMPLEX, INTENT(IN)     :: cpw2(ngpt2)
      REAL, INTENT(IN)        :: olappw(ngpt1, ngpt2)
      REAL, INTENT(IN)        :: olapmt(hybrid%maxindx, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype)
Daniel Wortmann's avatar
Daniel Wortmann committed
219
!     - local -
Matthias Redies's avatar
Matthias Redies committed
220
221
222
      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
223
224

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

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

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

!       wfolap = wfolap + cdum

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

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

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

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

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

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

!       wfolap = wfolap + cdum

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

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

Matthias Redies's avatar
Matthias Redies committed
313
      USE m_types
314
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
315
316
317
318
319

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

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

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

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

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

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

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

Matthias Redies's avatar
Matthias Redies committed
368
369
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms
370

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

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

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

   END FUNCTION wfolap2
END MODULE m_olap