olap.F90 10.9 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1
MODULE m_olap
2
      USE m_types_hybdat
3

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

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

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

Matthias Redies's avatar
Matthias Redies committed
29
      call timestart("olap_pw")
Matthias Redies's avatar
Matthias Redies committed
30 31 32 33
      DO i = 1, ngpt
         DO j = 1, i
            dg = gpt(:, j) - gpt(:, i)
            g = gptnorm(dg, cell%bmat)
Matthias Redies's avatar
Matthias Redies committed
34
            IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
               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
Matthias Redies's avatar
Matthias Redies committed
51
                        olap%data_r(i, j) = real(olap%data_r(i, j) - fgr*exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent))))
Matthias Redies's avatar
Matthias Redies committed
52 53 54 55 56 57 58 59 60 61 62 63 64 65
                     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
66
      END DO
Matthias Redies's avatar
Matthias Redies committed
67
      call timestop("olap_pw")
Matthias Redies's avatar
Matthias Redies committed
68
   END SUBROUTINE olap_pw
69

Matthias Redies's avatar
Matthias Redies committed
70
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71 72 73

!     olap_pwp  calculates upper triangular part of overlap matrix

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

Matthias Redies's avatar
Matthias Redies committed
76
      USE m_constants, ONLY: REAL_NOT_INITALIZED, CMPLX_NOT_INITALIZED
77 78
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
79 80
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms
81 82

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

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

Matthias Redies's avatar
Matthias Redies committed
96 97
      olap_r=REAL_NOT_INITALIZED; olap_c=CMPLX_NOT_INITALIZED

Daniel Wortmann's avatar
Daniel Wortmann committed
98
      if (l_real) THEN
Matthias Redies's avatar
Matthias Redies committed
99 100 101 102 103 104 105
         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
106
               IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
107 108 109 110 111 112 113 114 115 116 117
                  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
118
                        olap_r(k) = olap_r(k) - real(fgr* &
119
                             exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent))))
Matthias Redies's avatar
Matthias Redies committed
120 121 122 123
                     END DO
                  END DO
               END IF
               IF (i == j) olap_r(k) = olap_r(k) + 1
124
            END DO
Matthias Redies's avatar
Matthias Redies committed
125 126 127 128 129 130 131 132 133
         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
134
               IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
135 136 137 138 139 140 141 142 143 144 145 146
                  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* &
147
                             exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
Matthias Redies's avatar
Matthias Redies committed
148 149 150 151
                     END DO
                  END DO
               END IF
               IF (i == j) olap_c(k) = olap_c(k) + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
152
            END DO
Matthias Redies's avatar
Matthias Redies committed
153
         END DO
154

Matthias Redies's avatar
Matthias Redies committed
155 156
      endif
   END SUBROUTINE olap_pwp
157 158 159

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

Matthias Redies's avatar
Matthias Redies committed
160
   SUBROUTINE wfolap_init(olappw, olapmt, gpt,&
161
                          atoms, mpdata, cell,&
162
                          bas1, bas2)
Matthias Redies's avatar
Matthias Redies committed
163

Matthias Redies's avatar
Matthias Redies committed
164
      USE m_intgrf, ONLY: intgrf, intgrf_init
165 166
      USE m_types
      IMPLICIT NONE
167
      TYPE(t_mpdata), intent(in) :: mpdata
Matthias Redies's avatar
Matthias Redies committed
168 169 170
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms

171
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
172
      INTEGER, INTENT(IN)       :: gpt(:, :)!(3,ngpt)
173 174
      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)
Matthias Redies's avatar
Matthias Redies committed
175 176 177 178
      REAL, INTENT(INOUT)      :: olapmt(maxval(mpdata%num_radfun_per_l), &
                                          maxval(mpdata%num_radfun_per_l), &
                                          0:atoms%lmaxd, &
                                          atoms%ntype)
Matthias Redies's avatar
Matthias Redies committed
179
      TYPE(t_mat), INTENT(INOUT):: olappw
180 181

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

Matthias Redies's avatar
Matthias Redies committed
184
      REAL, ALLOCATABLE         :: gridf(:, :)
185

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

Matthias Redies's avatar
Matthias Redies committed
204 205 206
      CALL olap_pw(olappw, gpt, size(gpt, 2), atoms, cell)

   END SUBROUTINE wfolap_init
207

208
   FUNCTION wfolap_inv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpdata)
209 210 211 212

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

Matthias Redies's avatar
Matthias Redies committed
216
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
217
      COMPLEX                :: wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
218
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
219 220 221 222
      COMPLEX, INTENT(IN)     :: cmt1(:,:), cmt2(:,:)
      REAL, INTENT(IN)        :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      REAL, INTENT(IN)        :: olappw(:,:)
223
      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
224
!     - local -
Matthias Redies's avatar
Matthias Redies committed
225
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
Daniel Wortmann's avatar
Daniel Wortmann committed
226 227

      wfolap_inv = 0
Matthias Redies's avatar
Matthias Redies committed
228 229 230 231 232 233 234
      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
235
                  nn = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
236
                  wfolap_inv = wfolap_inv + &
237 238 239
                           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
240 241
                  lm = lm + nn
               END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
242
            END DO
Matthias Redies's avatar
Matthias Redies committed
243
         END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
244
      END DO
Matthias Redies's avatar
Matthias Redies committed
245 246

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

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

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
257
   END FUNCTION wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
258

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

Matthias Redies's avatar
Matthias Redies committed
261

Daniel Wortmann's avatar
Daniel Wortmann committed
262 263 264
      USE m_wrapper
      USE m_types
      IMPLICIT NONE
265
      TYPE(t_mpdata), intent(in) :: mpdata
Matthias Redies's avatar
Matthias Redies committed
266
      TYPE(t_atoms), INTENT(IN)   :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
267

Matthias Redies's avatar
Matthias Redies committed
268
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
269
      COMPLEX                :: wfolap_noinv
Matthias Redies's avatar
Matthias Redies committed
270
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
271 272 273 274
      COMPLEX, INTENT(IN)     :: cmt1(:,:),cmt2(:,:)
      COMPLEX, INTENT(IN)     :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      COMPLEX, INTENT(IN)     :: olappw(:,:)
275
      REAL, INTENT(IN)        :: olapmt(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
276
!     - local -
Matthias Redies's avatar
Matthias Redies committed
277
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
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 = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
288
                  wfolap_noinv = wfolap_noinv + &
289 290 291
                           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
292 293
                  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 310
   END FUNCTION wfolap_noinv
END MODULE m_olap