olap.F90 10.9 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

      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
32
            IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
               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
49
                        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
50 51 52 53 54 55 56 57 58 59 60 61 62 63
                     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
         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
102
               IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
103 104 105 106 107 108 109 110 111 112 113
                  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
114 115
                        olap_r(k) = olap_r(k) - real(fgr* &
             &               exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent))))
Matthias Redies's avatar
Matthias Redies committed
116 117 118 119
                     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
         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
130
               IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
                  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
   SUBROUTINE wfolap_init(olappw, olapmt, gpt,&
157
  &                       atoms, mpbasis, hybrid, cell,&
Matthias Redies's avatar
Matthias Redies committed
158 159
  &                       bas1, bas2)

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

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

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

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

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

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

   END SUBROUTINE wfolap_init
201

202
   FUNCTION wfolap_inv(cmt1, cpw1, cmt2, cpw2, ngpt1, ngpt2, olappw, olapmt, atoms, mpbasis,hybrid)
203 204 205 206

      USE m_wrapper
      USE m_types
      IMPLICIT NONE
207
      TYPE(t_mpbasis), intent(in) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
208 209
      TYPE(t_hybrid), INTENT(IN)   :: hybrid
      TYPE(t_atoms), INTENT(IN)   :: atoms
210

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

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

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

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

!       wfolap = wfolap + cdum

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

   FUNCTION wfolap_noinv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpbasis)

Matthias Redies's avatar
Matthias Redies committed
259

Daniel Wortmann's avatar
Daniel Wortmann committed
260 261 262
      USE m_wrapper
      USE m_types
      IMPLICIT NONE
263
      TYPE(t_mpbasis), intent(in) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
264
      TYPE(t_atoms), INTENT(IN)   :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
265

Matthias Redies's avatar
Matthias Redies committed
266
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
267
      COMPLEX                :: wfolap_noinv
Matthias Redies's avatar
Matthias Redies committed
268
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
269 270 271 272
      COMPLEX, INTENT(IN)     :: cmt1(:,:),cmt2(:,:)
      COMPLEX, INTENT(IN)     :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      COMPLEX, INTENT(IN)     :: olappw(:,:)
273
      REAL, INTENT(IN)        :: olapmt(maxval(mpbasis%num_radfun_per_l), maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
274
!     - local -
Matthias Redies's avatar
Matthias Redies committed
275
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
Matthias Redies's avatar
Matthias Redies committed
276
      REAL                   :: rdum1, rdum2
Daniel Wortmann's avatar
Daniel Wortmann committed
277
      wfolap_noinv = 0
Matthias Redies's avatar
Matthias Redies committed
278 279 280 281 282 283 284
      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
285
                  nn = mpbasis%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
286 287 288 289 290 291
                  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
292
            END DO
Matthias Redies's avatar
Matthias Redies committed
293
         END DO
294
      END DO
Matthias Redies's avatar
Matthias Redies committed
295 296

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

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

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
307 308
   END FUNCTION wfolap_noinv
END MODULE m_olap