olap.F90 10.6 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
                        olap_r(k) = olap_r(k) - real(fgr* &
115
                             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
                  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* &
143
                             exp(img*tpi_const*dot_product(dg, atoms%taual(:, icent)))
Matthias Redies's avatar
Matthias Redies committed
144 145 146 147
                     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 158
                          atoms, mpbasis, cell,&
                          bas1, bas2)
Matthias Redies's avatar
Matthias Redies committed
159

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
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms

167
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
168
      INTEGER, INTENT(IN)       :: gpt(:, :)!(3,ngpt)
169
      REAL, INTENT(IN)         ::  bas1(atoms%jmtd, maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype),&
170
                                  bas2(atoms%jmtd, maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
171
      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
172
      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
      DO itype = 1, atoms%ntype
         DO l = 0, atoms%lmax(itype)
183
            nn = mpbasis%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
184 185 186 187
            DO n2 = 1, nn
               DO n1 = 1, nn!n2
                  !IF( n1 .gt. 2 .or. n2 .gt. 2) CYCLE
                  olapmt(n1, n2, l, itype) = intgrf( &
188 189 190
                                      bas1(:, n1, l, itype)*bas1(:, n2, l, itype)&
                                     + bas2(:, n1, l, itype)*bas2(:, n2, l, itype),&
                                      atoms, 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

201
   FUNCTION wfolap_inv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpbasis)
202 203 204 205

      USE m_wrapper
      USE m_types
      IMPLICIT NONE
206
      TYPE(t_mpbasis), intent(in) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
207
      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
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
212 213 214 215
      COMPLEX, INTENT(IN)     :: cmt1(:,:), cmt2(:,:)
      REAL, INTENT(IN)        :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      REAL, INTENT(IN)        :: olappw(:,:)
216
      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
217
!     - local -
Matthias Redies's avatar
Matthias Redies committed
218
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
Daniel Wortmann's avatar
Daniel Wortmann committed
219 220

      wfolap_inv = 0
Matthias Redies's avatar
Matthias Redies committed
221 222 223 224 225 226 227
      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
228
                  nn = mpbasis%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
229
                  wfolap_inv = wfolap_inv + &
230 231 232
                           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
233 234
                  lm = lm + nn
               END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
235
            END DO
Matthias Redies's avatar
Matthias Redies committed
236
         END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
237
      END DO
Matthias Redies's avatar
Matthias Redies committed
238 239

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

241 242
!       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
243
!
244 245
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
Daniel Wortmann's avatar
Daniel Wortmann committed
246 247 248 249
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
250
   END FUNCTION wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
251 252 253

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

Matthias Redies's avatar
Matthias Redies committed
254

Daniel Wortmann's avatar
Daniel Wortmann committed
255 256 257
      USE m_wrapper
      USE m_types
      IMPLICIT NONE
258
      TYPE(t_mpbasis), intent(in) :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
259
      TYPE(t_atoms), INTENT(IN)   :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
260

Matthias Redies's avatar
Matthias Redies committed
261
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
262
      COMPLEX                :: wfolap_noinv
Matthias Redies's avatar
Matthias Redies committed
263
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
264 265 266 267
      COMPLEX, INTENT(IN)     :: cmt1(:,:),cmt2(:,:)
      COMPLEX, INTENT(IN)     :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      COMPLEX, INTENT(IN)     :: olappw(:,:)
268
      REAL, INTENT(IN)        :: olapmt(maxval(mpbasis%num_radfun_per_l), maxval(mpbasis%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
269
!     - local -
Matthias Redies's avatar
Matthias Redies committed
270
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
271

Daniel Wortmann's avatar
Daniel Wortmann committed
272
      wfolap_noinv = 0
Matthias Redies's avatar
Matthias Redies committed
273 274 275 276 277 278 279
      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
280
                  nn = mpbasis%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
281
                  wfolap_noinv = wfolap_noinv + &
282 283 284
                           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
285 286
                  lm = lm + nn
               END DO
287
            END DO
Matthias Redies's avatar
Matthias Redies committed
288
         END DO
289
      END DO
Matthias Redies's avatar
Matthias Redies committed
290 291

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

293 294
!       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
295
!
296 297
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
298 299 300 301
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
302 303
   END FUNCTION wfolap_noinv
END MODULE m_olap