olap.F90 11.1 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
      INTEGER                  :: i, j, itype, icent, ineq
      REAL                     :: g, r, fgr
Matthias Redies's avatar
Matthias Redies committed
26
      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 30 31 32 33 34
      if(olap%l_real) then
         olap%data_r = 0
      else 
         olap%data_c = 0
      endif
      
Matthias Redies's avatar
Matthias Redies committed
35
      call timestart("olap_pw")
Matthias Redies's avatar
Matthias Redies committed
36 37 38 39
      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
40
            IF (abs(g) < 1e-10) THEN
Matthias Redies's avatar
Matthias Redies committed
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
               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
57
                        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
58 59 60 61 62 63 64 65 66 67 68 69 70 71
                     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
72
      END DO
Matthias Redies's avatar
Matthias Redies committed
73
      call timestop("olap_pw")
Matthias Redies's avatar
Matthias Redies committed
74
   END SUBROUTINE olap_pw
75

Matthias Redies's avatar
Matthias Redies committed
76
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77 78 79

!     olap_pwp  calculates upper triangular part of overlap matrix

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

Matthias Redies's avatar
Matthias Redies committed
82 83
      USE m_constants, ONLY: REAL_NOT_INITALIZED, CMPLX_NOT_INITALIZED,&
                             fpi_const, tpi_const
84 85
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
86 87
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms
88 89

!     - scalars -
Matthias Redies's avatar
Matthias Redies committed
90
      INTEGER, INTENT(IN)       :: ngpt
91
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
92
      INTEGER, INTENT(IN)       :: gpt(:,:)
93

Matthias Redies's avatar
Matthias Redies committed
94
      LOGICAL, INTENT(IN)       :: l_real
Matthias Redies's avatar
Matthias Redies committed
95 96
      REAL, INTENT(INOUT)       ::  olap_r(ngpt*(ngpt + 1)/2)
      COMPLEX, INTENT(INOUT)    ::  olap_c(ngpt*(ngpt + 1)/2)
97
!     - local -
Matthias Redies's avatar
Matthias Redies committed
98 99 100
      INTEGER                  :: i, j, k, itype, icent, ineq
      REAL                     :: g, r, fgr
      COMPLEX, PARAMETER        :: img = (0.0, 1.0)
101 102
      INTEGER                  :: dg(3)

Matthias Redies's avatar
Matthias Redies committed
103 104
      olap_r=REAL_NOT_INITALIZED; olap_c=CMPLX_NOT_INITALIZED

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

Matthias Redies's avatar
Matthias Redies committed
162 163
      endif
   END SUBROUTINE olap_pwp
164 165 166

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

Matthias Redies's avatar
Matthias Redies committed
167
   SUBROUTINE wfolap_init(olappw, olapmt, gpt,&
168
                          atoms, mpdata, cell,&
169
                          bas1, bas2)
Matthias Redies's avatar
Matthias Redies committed
170

Matthias Redies's avatar
Matthias Redies committed
171
      USE m_intgrf, ONLY: intgrf, intgrf_init
172 173
      USE m_types
      IMPLICIT NONE
174
      TYPE(t_mpdata), intent(in) :: mpdata
Matthias Redies's avatar
Matthias Redies committed
175 176 177
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms

178
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
179
      INTEGER, INTENT(IN)       :: gpt(:, :)!(3,ngpt)
180 181
      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
182 183 184 185
      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
186
      TYPE(t_mat), INTENT(INOUT):: olappw
187 188

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

Matthias Redies's avatar
Matthias Redies committed
191
      REAL, ALLOCATABLE         :: gridf(:, :)
192

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

Matthias Redies's avatar
Matthias Redies committed
211 212 213
      CALL olap_pw(olappw, gpt, size(gpt, 2), atoms, cell)

   END SUBROUTINE wfolap_init
214

215
   FUNCTION wfolap_inv(cmt1, cpw1, cmt2, cpw2, olappw, olapmt, atoms, mpdata)
216 217 218 219

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

Matthias Redies's avatar
Matthias Redies committed
223
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
224
      COMPLEX                :: wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
225
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
226 227 228 229
      COMPLEX, INTENT(IN)     :: cmt1(:,:), cmt2(:,:)
      REAL, INTENT(IN)        :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      REAL, INTENT(IN)        :: olappw(:,:)
230
      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
231
!     - local -
Matthias Redies's avatar
Matthias Redies committed
232
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
Daniel Wortmann's avatar
Daniel Wortmann committed
233 234

      wfolap_inv = 0
Matthias Redies's avatar
Matthias Redies committed
235 236 237 238 239 240 241
      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
242
                  nn = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
243
                  wfolap_inv = wfolap_inv + &
244 245 246
                           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
247 248
                  lm = lm + nn
               END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
249
            END DO
Matthias Redies's avatar
Matthias Redies committed
250
         END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
251
      END DO
Matthias Redies's avatar
Matthias Redies committed
252 253

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

255 256
!       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
257
!
258 259
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
Daniel Wortmann's avatar
Daniel Wortmann committed
260 261 262 263
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
264
   END FUNCTION wfolap_inv
Matthias Redies's avatar
Matthias Redies committed
265

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

Matthias Redies's avatar
Matthias Redies committed
268

Daniel Wortmann's avatar
Daniel Wortmann committed
269 270 271
      USE m_wrapper
      USE m_types
      IMPLICIT NONE
272
      TYPE(t_mpdata), intent(in) :: mpdata
Matthias Redies's avatar
Matthias Redies committed
273
      TYPE(t_atoms), INTENT(IN)   :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
274

Matthias Redies's avatar
Matthias Redies committed
275
!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
276
      COMPLEX                :: wfolap_noinv
Matthias Redies's avatar
Matthias Redies committed
277
!     - arrays -
Matthias Redies's avatar
Matthias Redies committed
278 279 280 281
      COMPLEX, INTENT(IN)     :: cmt1(:,:),cmt2(:,:)
      COMPLEX, INTENT(IN)     :: cpw1(:)
      COMPLEX, INTENT(IN)     :: cpw2(:)
      COMPLEX, INTENT(IN)     :: olappw(:,:)
282
      REAL, INTENT(IN)        :: olapmt(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype)
283
!     - local -
Matthias Redies's avatar
Matthias Redies committed
284
      INTEGER                :: itype, ieq, iatom, l, m, lm, nn
285

Daniel Wortmann's avatar
Daniel Wortmann committed
286
      wfolap_noinv = 0
Matthias Redies's avatar
Matthias Redies committed
287 288 289 290 291 292 293
      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
294
                  nn = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
295
                  wfolap_noinv = wfolap_noinv + &
296 297 298
                           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
299 300
                  lm = lm + nn
               END DO
301
            END DO
Matthias Redies's avatar
Matthias Redies committed
302
         END DO
303
      END DO
Matthias Redies's avatar
Matthias Redies committed
304 305

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

307 308
!       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
309
!
310 311
!       rdum1 = dot_product(cpw1,rarr1)
!       rdum2 = dot_product(cpw1,rarr2)
312 313 314 315
!       cdum  = cmplx( rdum1, rdum2 )

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
316 317
   END FUNCTION wfolap_noinv
END MODULE m_olap