olap.F90 13.6 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1
MODULE m_olap
2
   USE m_types_hybdat
Matthias Redies's avatar
Matthias Redies committed
3
#ifdef CPP_MPI
4
   use mpi
Matthias Redies's avatar
Matthias Redies committed
5
#endif
6 7
   private olap_pw_real, olap_pw_cmplx
   public olap_pw, olap_pwp, wfolap_init, wfolap_inv, wfolap_noinv
8

Matthias Redies's avatar
Matthias Redies committed
9
CONTAINS
10 11 12 13 14

!     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

15
   SUBROUTINE olap_pw(olap, gpt, ngpt, atoms, cell, fmpi)
Matthias Redies's avatar
Matthias Redies committed
16
      use m_juDFT
Matthias Redies's avatar
Matthias Redies committed
17
      USE m_constants
18 19
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
20 21
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms
22
      type(t_mpi), intent(in)    :: fmpi
23 24

!     - scalars -
Matthias Redies's avatar
Matthias Redies committed
25
      INTEGER, INTENT(IN)       :: ngpt
26
!     - arrays -
27
      INTEGER, INTENT(IN)       :: gpt(:, :)
Daniel Wortmann's avatar
Daniel Wortmann committed
28
      TYPE(t_mat)              :: olap
Matthias Redies's avatar
Matthias Redies committed
29

30
      if (olap%l_real) then
31
         call olap_pw_real(olap, gpt, ngpt, atoms, cell, fmpi)
32 33
      else
         call olap_pw_cmplx(olap, gpt, ngpt, atoms, cell)
Matthias Redies's avatar
Matthias Redies committed
34
      endif
35 36
   END SUBROUTINE olap_pw

37
   subroutine olap_pw_real(olap, gpt, ngpt, atoms, cell, fmpi)
38 39 40 41 42
      use m_juDFT
      USE m_constants
      USE m_types
      IMPLICIT NONE
      TYPE(t_cell), INTENT(IN)   :: cell
43 44
      TYPE(t_atoms), INTENT(IN)  :: atoms
      type(t_mpi), intent(in)    :: fmpi
45 46 47 48 49 50 51

      !     - scalars -
      INTEGER, INTENT(IN)       :: ngpt
      !     - arrays -
      INTEGER, INTENT(IN)       :: gpt(:, :)
      TYPE(t_mat)              :: olap
      !     - local -
52
      INTEGER                  :: i, j, itype, icent, ierr, root
53 54 55 56 57
      REAL                     :: g, r, fgr
      COMPLEX, PARAMETER       :: img = (0.0, 1.0)
      INTEGER                  :: dg(3)

      call timestart("olap_pw_real")
58
      !$OMP PARALLEL default(none) &
59
      !$OMP private(i,j,dg,g,itype, icent, r, fgr)&
60
      !$OMP shared(olap, gpt, cell, atoms, ngpt, fmpi)
61 62

      !$OMP DO schedule(guided) 
63
      DO j = 1+fmpi%n_rank, ngpt, fmpi%n_size
64
         DO i = 1, j
65
            olap%data_r(i,j) = 0.0
Matthias Redies's avatar
Matthias Redies committed
66 67
            dg = gpt(:, j) - gpt(:, i)
            g = gptnorm(dg, cell%bmat)
Matthias Redies's avatar
Matthias Redies committed
68
            IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
69 70
               DO itype = 1, atoms%ntype
                  r = atoms%rmt(itype)
71
                  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
72 73
               END DO
            ELSE
74 75 76 77 78 79 80 81 82
               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
         END DO
      END DO
83 84 85 86
      !$OMP end do

      ! work on diagonal
      !$OMP DO
87
      do i = 1+fmpi%n_rank, ngpt, fmpi%n_size
88 89 90 91
         olap%data_r(i,i) = olap%data_r(i,i) + 1
      enddo
      !$OMP END DO
      !$OMP end parallel
92 93 94 95 96 97
#ifdef CPP_MPI
      do j = 1, ngpt
         root = mod(j-1, fmpi%n_size)
         call MPI_Bcast(olap%data_r(1,j), j, MPI_DOUBLE_PRECISION, root, fmpi%sub_comm, ierr)
      enddo
#endif
98
      call olap%u2l()
99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
      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 -
116
      INTEGER                  :: i, j, itype, icent
117 118 119 120 121 122 123
      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
124
      !$OMP shared(olap, gpt, cell, atoms, ngpt)
125 126 127 128 129 130
      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
131
               DO itype = 1, atoms%ntype
132 133 134 135 136 137 138
                  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
139 140
                  r = g*atoms%rmt(itype)
                  fgr = fpi_const*(sin(r) - r*cos(r))/g**3/cell%omtil
141
                  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
142 143
               END DO
            END IF
144 145
            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
146
         END DO
147
      END DO
148 149 150
      !$OMP END PARALLEL DO
      call timestop("olap_pw_cmplx")
   END SUBROUTINE olap_pw_cmplx
151

Matthias Redies's avatar
Matthias Redies committed
152
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153 154 155

!     olap_pwp  calculates upper triangular part of overlap matrix

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

158
      USE m_constants, ONLY: REAL_NOT_INITALIZED, CMPLX_NOT_INITALIZED, &
Matthias Redies's avatar
Matthias Redies committed
159
                             fpi_const, tpi_const
160 161
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
162 163
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms
164 165

!     - scalars -
Matthias Redies's avatar
Matthias Redies committed
166
      INTEGER, INTENT(IN)       :: ngpt
167
!     - arrays -
168
      INTEGER, INTENT(IN)       :: gpt(:, :)
169

Matthias Redies's avatar
Matthias Redies committed
170
      LOGICAL, INTENT(IN)       :: l_real
Matthias Redies's avatar
Matthias Redies committed
171 172
      REAL, INTENT(INOUT)       ::  olap_r(ngpt*(ngpt + 1)/2)
      COMPLEX, INTENT(INOUT)    ::  olap_c(ngpt*(ngpt + 1)/2)
173
!     - local -
Matthias Redies's avatar
Matthias Redies committed
174 175 176
      INTEGER                  :: i, j, k, itype, icent, ineq
      REAL                     :: g, r, fgr
      COMPLEX, PARAMETER        :: img = (0.0, 1.0)
177 178
      INTEGER                  :: dg(3)

179
      olap_r = REAL_NOT_INITALIZED; olap_c = CMPLX_NOT_INITALIZED
Matthias Redies's avatar
Matthias Redies committed
180

Daniel Wortmann's avatar
Daniel Wortmann committed
181
      if (l_real) THEN
Matthias Redies's avatar
Matthias Redies committed
182 183 184 185 186 187 188
         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
189
               IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
190 191 192 193 194 195 196 197 198 199 200
                  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
201
                        olap_r(k) = olap_r(k) - real(fgr* &
202
                                                     exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent))))
Matthias Redies's avatar
Matthias Redies committed
203 204 205 206
                     END DO
                  END DO
               END IF
               IF (i == j) olap_r(k) = olap_r(k) + 1
207
            END DO
Matthias Redies's avatar
Matthias Redies committed
208 209 210 211 212 213 214 215 216
         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
217
               IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
218 219 220 221 222 223 224 225 226 227 228 229
                  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* &
230
                                    exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
Matthias Redies's avatar
Matthias Redies committed
231 232 233 234
                     END DO
                  END DO
               END IF
               IF (i == j) olap_c(k) = olap_c(k) + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
235
            END DO
Matthias Redies's avatar
Matthias Redies committed
236
         END DO
237

Matthias Redies's avatar
Matthias Redies committed
238 239
      endif
   END SUBROUTINE olap_pwp
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
!    SUBROUTINE wfolap_init(olappw, olapmt, gpt, &
!                           atoms, mpdata, cell, &
!                           bas1, bas2)

!       USE m_intgrf, ONLY: intgrf, intgrf_init
!       USE m_types
!       IMPLICIT NONE
!       TYPE(t_mpdata), intent(in) :: mpdata
!       TYPE(t_cell), INTENT(IN)   :: cell
!       TYPE(t_atoms), INTENT(IN)   :: atoms

! !     - arrays -
!       INTEGER, INTENT(IN)       :: gpt(:, :)!(3,ngpt)
!       REAL, INTENT(IN)         ::  bas1(atoms%jmtd, maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype), &
!                                   bas2(atoms%jmtd, maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
!       REAL, INTENT(INOUT)      :: olapmt(maxval(mpdata%num_radfun_per_l), &
!                                          maxval(mpdata%num_radfun_per_l), &
!                                          0:atoms%lmaxd, &
!                                          atoms%ntype)
!       TYPE(t_mat), INTENT(INOUT):: olappw

! !     - local -
!       INTEGER                  :: itype, l, nn, n1, n2

!       REAL, ALLOCATABLE         :: gridf(:, :)

!       CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
!       olapmt = 0
!       DO itype = 1, atoms%ntype
!          DO l = 0, atoms%lmax(itype)
!             nn = mpdata%num_radfun_per_l(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, itype, gridf)
! !               olapmt(n2,n1,l,itype) = olapmt(n1,n2,l,itype)
!                END DO
!             END DO
!          END DO
!       END DO

!       CALL olap_pw(olappw, gpt, size(gpt, 2), atoms, cell)

!    END SUBROUTINE wfolap_init
290

291
   FUNCTION wfolap_inv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpdata)
292 293 294 295

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

Matthias Redies's avatar
Matthias Redies committed
299
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
300
      COMPLEX                :: wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
301
!     - arrays -
302
      COMPLEX, INTENT(IN)     :: cmt1(:, :), cmt2(:, :)
Matthias Redies's avatar
Matthias Redies committed
303 304
      REAL, INTENT(IN)        :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
305
      REAL, INTENT(IN)        :: olappw(:, :)
306
      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
307
!     - local -
Matthias Redies's avatar
Matthias Redies committed
308
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
Daniel Wortmann's avatar
Daniel Wortmann committed
309 310

      wfolap_inv = 0
Matthias Redies's avatar
Matthias Redies committed
311 312 313 314 315 316 317
      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
318
                  nn = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
319
                  wfolap_inv = wfolap_inv + &
320 321 322
                               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
323 324
                  lm = lm + nn
               END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
325
            END DO
Matthias Redies's avatar
Matthias Redies committed
326
         END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
327
      END DO
Matthias Redies's avatar
Matthias Redies committed
328 329

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

331 332
!       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
333
!
334 335
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
Daniel Wortmann's avatar
Daniel Wortmann committed
336 337 338 339
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
340
   END FUNCTION wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
341

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

Daniel Wortmann's avatar
Daniel Wortmann committed
344 345 346
      USE m_wrapper
      USE m_types
      IMPLICIT NONE
347
      TYPE(t_mpdata), intent(in) :: mpdata
Matthias Redies's avatar
Matthias Redies committed
348
      TYPE(t_atoms), INTENT(IN)   :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
349

Matthias Redies's avatar
Matthias Redies committed
350
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
351
      COMPLEX                :: wfolap_noinv
Matthias Redies's avatar
Matthias Redies committed
352
!     - arrays -
353
      COMPLEX, INTENT(IN)     :: cmt1(:, :), cmt2(:, :)
Matthias Redies's avatar
Matthias Redies committed
354 355
      COMPLEX, INTENT(IN)     :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
356
      COMPLEX, INTENT(IN)     :: olappw(:, :)
357
      REAL, INTENT(IN)        :: olapmt(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
358
!     - local -
Matthias Redies's avatar
Matthias Redies committed
359
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
360

Daniel Wortmann's avatar
Daniel Wortmann committed
361
      wfolap_noinv = 0
Matthias Redies's avatar
Matthias Redies committed
362 363 364 365 366 367 368
      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
369
                  nn = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
370
                  wfolap_noinv = wfolap_noinv + &
371 372 373
                                 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
374 375
                  lm = lm + nn
               END DO
376
            END DO
Matthias Redies's avatar
Matthias Redies committed
377
         END DO
378
      END DO
Matthias Redies's avatar
Matthias Redies committed
379 380

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

382 383
!       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
384
!
385 386
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
387 388 389 390
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
391 392
   END FUNCTION wfolap_noinv
END MODULE m_olap