olap.F90 10.8 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
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 88 89
      LOGICAL, INTENT(IN)       :: l_real
      REAL, INTENT(OUT)         ::  olap_r(ngpt*(ngpt + 1)/2)
      COMPLEX, INTENT(OUT)      ::  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)

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

Matthias Redies's avatar
Matthias Redies committed
153 154
      endif
   END SUBROUTINE olap_pwp
155 156 157

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

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

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

169
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
170
      INTEGER, INTENT(IN)       :: gpt(:, :)!(3,ngpt)
171 172
      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
173 174 175 176
      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
177
      TYPE(t_mat), INTENT(INOUT):: olappw
178 179

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

Matthias Redies's avatar
Matthias Redies committed
182
      REAL, ALLOCATABLE         :: gridf(:, :)
183

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

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

   END SUBROUTINE wfolap_init
205

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

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

Matthias Redies's avatar
Matthias Redies committed
214
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
215
      COMPLEX                :: wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
216
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
217 218 219 220
      COMPLEX, INTENT(IN)     :: cmt1(:,:), cmt2(:,:)
      REAL, INTENT(IN)        :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      REAL, INTENT(IN)        :: olappw(:,:)
221
      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
222
!     - local -
Matthias Redies's avatar
Matthias Redies committed
223
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
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 = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
234
                  wfolap_inv = wfolap_inv + &
235 236 237
                           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
238 239
                  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
   FUNCTION wfolap_noinv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpdata)
Matthias Redies's avatar
Matthias Redies committed
258

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_mpdata), intent(in) :: mpdata
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(mpdata%num_radfun_per_l), maxval(mpdata%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
276

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 = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
286
                  wfolap_noinv = wfolap_noinv + &
287 288 289
                           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
290 291
                  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