olap.F90 13.3 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1
MODULE m_olap
2 3 4
   USE m_types_hybdat
   private olap_pw_real, olap_pw_cmplx
   public olap_pw, olap_pwp, wfolap_init, wfolap_inv, wfolap_noinv
5

Matthias Redies's avatar
Matthias Redies committed
6
CONTAINS
7 8 9 10 11

!     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
12
   SUBROUTINE olap_pw(olap, gpt, ngpt, atoms, cell)
Matthias Redies's avatar
Matthias Redies committed
13
      use m_juDFT
Matthias Redies's avatar
Matthias Redies committed
14
      USE m_constants
15 16
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
17 18
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms
19 20

!     - scalars -
Matthias Redies's avatar
Matthias Redies committed
21
      INTEGER, INTENT(IN)       :: ngpt
22
!     - arrays -
23
      INTEGER, INTENT(IN)       :: gpt(:, :)
Daniel Wortmann's avatar
Daniel Wortmann committed
24
      TYPE(t_mat)              :: olap
25
!     - local -
Matthias Redies's avatar
Matthias Redies committed
26 27
      INTEGER                  :: i, j, itype, icent, ineq
      REAL                     :: g, r, fgr
Matthias Redies's avatar
Matthias Redies committed
28
      COMPLEX, PARAMETER       :: img = (0.0, 1.0)
29
      INTEGER                  :: dg(3)
Matthias Redies's avatar
Matthias Redies committed
30

31 32 33 34
      if (olap%l_real) then
         call olap_pw_real(olap, gpt, ngpt, atoms, cell)
      else
         call olap_pw_cmplx(olap, gpt, ngpt, atoms, cell)
Matthias Redies's avatar
Matthias Redies committed
35
      endif
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
   END SUBROUTINE olap_pw

   subroutine olap_pw_real(olap, gpt, ngpt, atoms, cell)
      use m_juDFT
      USE m_constants
      USE m_types
      IMPLICIT NONE
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms

      !     - scalars -
      INTEGER, INTENT(IN)       :: ngpt
      !     - arrays -
      INTEGER, INTENT(IN)       :: gpt(:, :)
      TYPE(t_mat)              :: olap
      !     - local -
      INTEGER                  :: i, j, itype, icent, ineq
      REAL                     :: g, r, fgr
      COMPLEX, PARAMETER       :: img = (0.0, 1.0)
      INTEGER                  :: dg(3)

      call timestart("olap_pw_real")
      !$OMP PARALLEL DO default(none) schedule(guided) &
      !$OMP private(i,j,dg,g,itype, icent, r, fgr)&
Matthias Redies's avatar
Matthias Redies committed
60
      !$OMP shared(olap, gpt, cell, atoms, ngpt)
Matthias Redies's avatar
Matthias Redies committed
61 62
      DO i = 1, ngpt
         DO j = 1, i
63
            olap%data_r(i,j) = 0.0
Matthias Redies's avatar
Matthias Redies committed
64 65
            dg = gpt(:, j) - gpt(:, i)
            g = gptnorm(dg, cell%bmat)
Matthias Redies's avatar
Matthias Redies committed
66
            IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
67 68
               DO itype = 1, atoms%ntype
                  r = atoms%rmt(itype)
69
                  olap%data_r(i, j) = olap%data_r(i, j) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
Matthias Redies's avatar
Matthias Redies committed
70 71
               END DO
            ELSE
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
               do icent = 1, atoms%nat
                  itype = atoms%itype(icent)
                  r = g*atoms%rmt(itype)
                  fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
                  olap%data_r(i, j) = real(olap%data_r(i, j) - fgr*exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent))))
               END DO
            END IF
            IF (i == j) olap%data_r(i, j) = olap%data_r(i, j) + 1
            olap%data_r(j, i) = olap%data_r(i, j)
         END DO
      END DO
      !$OMP END PARALLEL DO
      call timestop("olap_pw_real")
   END SUBROUTINE olap_pw_real

   SUBROUTINE olap_pw_cmplx(olap, gpt, ngpt, atoms, cell)
      use m_juDFT
      USE m_constants
      USE m_types
      IMPLICIT NONE
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms

!     - scalars -
      INTEGER, INTENT(IN)       :: ngpt
!     - arrays -
      INTEGER, INTENT(IN)       :: gpt(:, :)
      TYPE(t_mat)              :: olap
!     - local -
      INTEGER                  :: i, j, itype, icent, ineq
      REAL                     :: g, r, fgr
      COMPLEX, PARAMETER       :: img = (0.0, 1.0)
      INTEGER                  :: dg(3)

      call timestart("olap_pw_cmplx")
      !$OMP PARALLEL DO default(none) schedule(guided) &
      !$OMP private(i,j,dg,g,itype, icent, r, fgr)&
Matthias Redies's avatar
Matthias Redies committed
109
      !$OMP shared(olap, gpt, cell, atoms, ngpt)
110 111 112 113 114 115
      DO i = 1, ngpt
         DO j = 1, i
            olap%data_c(i,j) = cmplx_0
            dg = gpt(:, j) - gpt(:, i)
            g = gptnorm(dg, cell%bmat)
            IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
116
               DO itype = 1, atoms%ntype
117 118 119 120 121 122 123
                  r = atoms%rmt(itype)
                  olap%data_c(i, j) = olap%data_c(i, j) - atoms%neq(itype)*fpi_const*r**3/3/cell%omtil
               END DO
            ELSE
               icent = 0
               do icent = 1, atoms%nat
                  itype = atoms%itype(icent)
Matthias Redies's avatar
Matthias Redies committed
124 125
                  r = g*atoms%rmt(itype)
                  fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
126
                  olap%data_c(i, j) = olap%data_c(i, j) - fgr*exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
Matthias Redies's avatar
Matthias Redies committed
127 128
               END DO
            END IF
129 130
            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))
Matthias Redies's avatar
Matthias Redies committed
131
         END DO
132
      END DO
133 134 135
      !$OMP END PARALLEL DO
      call timestop("olap_pw_cmplx")
   END SUBROUTINE olap_pw_cmplx
136

Matthias Redies's avatar
Matthias Redies committed
137
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138 139 140

!     olap_pwp  calculates upper triangular part of overlap matrix

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

143
      USE m_constants, ONLY: REAL_NOT_INITALIZED, CMPLX_NOT_INITALIZED, &
Matthias Redies's avatar
Matthias Redies committed
144
                             fpi_const, tpi_const
145 146
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
147 148
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms
149 150

!     - scalars -
Matthias Redies's avatar
Matthias Redies committed
151
      INTEGER, INTENT(IN)       :: ngpt
152
!     - arrays -
153
      INTEGER, INTENT(IN)       :: gpt(:, :)
154

Matthias Redies's avatar
Matthias Redies committed
155
      LOGICAL, INTENT(IN)       :: l_real
Matthias Redies's avatar
Matthias Redies committed
156 157
      REAL, INTENT(INOUT)       ::  olap_r(ngpt*(ngpt + 1)/2)
      COMPLEX, INTENT(INOUT)    ::  olap_c(ngpt*(ngpt + 1)/2)
158
!     - local -
Matthias Redies's avatar
Matthias Redies committed
159 160 161
      INTEGER                  :: i, j, k, itype, icent, ineq
      REAL                     :: g, r, fgr
      COMPLEX, PARAMETER        :: img = (0.0, 1.0)
162 163
      INTEGER                  :: dg(3)

164
      olap_r = REAL_NOT_INITALIZED; olap_c = CMPLX_NOT_INITALIZED
Matthias Redies's avatar
Matthias Redies committed
165

Daniel Wortmann's avatar
Daniel Wortmann committed
166
      if (l_real) THEN
Matthias Redies's avatar
Matthias Redies committed
167 168 169 170 171 172 173
         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
Matthias Redies's avatar
Matthias Redies committed
174
               IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
175 176 177 178 179 180 181 182 183 184 185
                  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
Matthias Redies's avatar
Matthias Redies committed
186
                        olap_r(k) = olap_r(k) - real(fgr* &
187
                                                     exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent))))
Matthias Redies's avatar
Matthias Redies committed
188 189 190 191
                     END DO
                  END DO
               END IF
               IF (i == j) olap_r(k) = olap_r(k) + 1
192
            END DO
Matthias Redies's avatar
Matthias Redies committed
193 194 195 196 197 198 199 200 201
         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
Matthias Redies's avatar
Matthias Redies committed
202
               IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
203 204 205 206 207 208 209 210 211 212 213 214
                  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* &
215
                                    exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
Matthias Redies's avatar
Matthias Redies committed
216 217 218 219
                     END DO
                  END DO
               END IF
               IF (i == j) olap_c(k) = olap_c(k) + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
220
            END DO
Matthias Redies's avatar
Matthias Redies committed
221
         END DO
222

Matthias Redies's avatar
Matthias Redies committed
223 224
      endif
   END SUBROUTINE olap_pwp
225 226 227

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

228 229
   SUBROUTINE wfolap_init(olappw, olapmt, gpt, &
                          atoms, mpdata, cell, &
230
                          bas1, bas2)
Matthias Redies's avatar
Matthias Redies committed
231

Matthias Redies's avatar
Matthias Redies committed
232
      USE m_intgrf, ONLY: intgrf, intgrf_init
233 234
      USE m_types
      IMPLICIT NONE
235
      TYPE(t_mpdata), intent(in) :: mpdata
Matthias Redies's avatar
Matthias Redies committed
236 237 238
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms

239
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
240
      INTEGER, INTENT(IN)       :: gpt(:, :)!(3,ngpt)
241
      REAL, INTENT(IN)         ::  bas1(atoms%jmtd, maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype), &
242
                                  bas2(atoms%jmtd, maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
Matthias Redies's avatar
Matthias Redies committed
243
      REAL, INTENT(INOUT)      :: olapmt(maxval(mpdata%num_radfun_per_l), &
244 245 246
                                         maxval(mpdata%num_radfun_per_l), &
                                         0:atoms%lmaxd, &
                                         atoms%ntype)
Matthias Redies's avatar
Matthias Redies committed
247
      TYPE(t_mat), INTENT(INOUT):: olappw
248 249

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

Matthias Redies's avatar
Matthias Redies committed
252
      REAL, ALLOCATABLE         :: gridf(:, :)
253

Matthias Redies's avatar
Matthias Redies committed
254
      CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
255
      olapmt = 0
Matthias Redies's avatar
Matthias Redies committed
256 257
      DO itype = 1, atoms%ntype
         DO l = 0, atoms%lmax(itype)
258
            nn = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
259 260 261 262
            DO n2 = 1, nn
               DO n1 = 1, nn!n2
                  !IF( n1 .gt. 2 .or. n2 .gt. 2) CYCLE
                  olapmt(n1, n2, l, itype) = intgrf( &
263 264 265
                                             bas1(:, n1, l, itype)*bas1(:, n2, l, itype) &
                                             + bas2(:, n1, l, itype)*bas2(:, n2, l, itype), &
                                             atoms, itype, gridf)
266
!               olapmt(n2,n1,l,itype) = olapmt(n1,n2,l,itype)
Matthias Redies's avatar
Matthias Redies committed
267
               END DO
268
            END DO
Matthias Redies's avatar
Matthias Redies committed
269
         END DO
270 271
      END DO

Matthias Redies's avatar
Matthias Redies committed
272 273 274
      CALL olap_pw(olappw, gpt, size(gpt, 2), atoms, cell)

   END SUBROUTINE wfolap_init
275

276
   FUNCTION wfolap_inv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpdata)
277 278 279 280

      USE m_wrapper
      USE m_types
      IMPLICIT NONE
281
      TYPE(t_mpdata), intent(in) :: mpdata
Matthias Redies's avatar
Matthias Redies committed
282
      TYPE(t_atoms), INTENT(IN)   :: atoms
283

Matthias Redies's avatar
Matthias Redies committed
284
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
285
      COMPLEX                :: wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
286
!     - arrays -
287
      COMPLEX, INTENT(IN)     :: cmt1(:, :), cmt2(:, :)
Matthias Redies's avatar
Matthias Redies committed
288 289
      REAL, INTENT(IN)        :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
290
      REAL, INTENT(IN)        :: olappw(:, :)
291
      REAL, INTENT(IN)        :: olapmt(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
Daniel Wortmann's avatar
Daniel Wortmann committed
292
!     - local -
Matthias Redies's avatar
Matthias Redies committed
293
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
Daniel Wortmann's avatar
Daniel Wortmann committed
294 295

      wfolap_inv = 0
Matthias Redies's avatar
Matthias Redies committed
296 297 298 299 300 301 302
      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
303
                  nn = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
304
                  wfolap_inv = wfolap_inv + &
305 306 307
                               dot_product(cmt1(lm + 1:lm + nn, iatom), &
                                           matmul(olapmt(:nn, :nn, l, itype), &
                                                  cmt2(lm + 1:lm + nn, iatom)))
Matthias Redies's avatar
Matthias Redies committed
308 309
                  lm = lm + nn
               END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
310
            END DO
Matthias Redies's avatar
Matthias Redies committed
311
         END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
312
      END DO
Matthias Redies's avatar
Matthias Redies committed
313 314

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

316 317
!       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
318
!
319 320
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
Daniel Wortmann's avatar
Daniel Wortmann committed
321 322 323 324
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
325
   END FUNCTION wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
326

327
   FUNCTION wfolap_noinv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpdata)
Matthias Redies's avatar
Matthias Redies committed
328

Daniel Wortmann's avatar
Daniel Wortmann committed
329 330 331
      USE m_wrapper
      USE m_types
      IMPLICIT NONE
332
      TYPE(t_mpdata), intent(in) :: mpdata
Matthias Redies's avatar
Matthias Redies committed
333
      TYPE(t_atoms), INTENT(IN)   :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
334

Matthias Redies's avatar
Matthias Redies committed
335
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
336
      COMPLEX                :: wfolap_noinv
Matthias Redies's avatar
Matthias Redies committed
337
!     - arrays -
338
      COMPLEX, INTENT(IN)     :: cmt1(:, :), cmt2(:, :)
Matthias Redies's avatar
Matthias Redies committed
339 340
      COMPLEX, INTENT(IN)     :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
341
      COMPLEX, INTENT(IN)     :: olappw(:, :)
342
      REAL, INTENT(IN)        :: olapmt(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
343
!     - local -
Matthias Redies's avatar
Matthias Redies committed
344
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
345

Daniel Wortmann's avatar
Daniel Wortmann committed
346
      wfolap_noinv = 0
Matthias Redies's avatar
Matthias Redies committed
347 348 349 350 351 352 353
      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
354
                  nn = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
355
                  wfolap_noinv = wfolap_noinv + &
356 357 358
                                 dot_product(cmt1(lm + 1:lm + nn, iatom), &
                                             matmul(olapmt(:nn, :nn, l, itype), &
                                                    cmt2(lm + 1:lm + nn, iatom)))
Matthias Redies's avatar
Matthias Redies committed
359 360
                  lm = lm + nn
               END DO
361
            END DO
Matthias Redies's avatar
Matthias Redies committed
362
         END DO
363
      END DO
Matthias Redies's avatar
Matthias Redies committed
364 365

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

367 368
!       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
369
!
370 371
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
372 373 374 375
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
376 377
   END FUNCTION wfolap_noinv
END MODULE m_olap