olap.F90 11 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 77
      USE m_constants, ONLY: REAL_NOT_INITALIZED, CMPLX_NOT_INITALIZED,&
                             fpi_const, tpi_const
78 79
      USE m_types
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
80 81
      TYPE(t_cell), INTENT(IN)   :: cell
      TYPE(t_atoms), INTENT(IN)   :: atoms
82 83

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

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

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

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

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

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

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

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

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

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

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

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

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

   END SUBROUTINE wfolap_init
208

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

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

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

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

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

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

!       wfolap = wfolap + cdum

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

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

Matthias Redies's avatar
Matthias Redies committed
262

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

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

Daniel Wortmann's avatar
Daniel Wortmann committed
280
      wfolap_noinv = 0
Matthias Redies's avatar
Matthias Redies committed
281 282 283 284 285 286 287
      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
288
                  nn = mpdata%num_radfun_per_l(l, itype)
Matthias Redies's avatar
Matthias Redies committed
289
                  wfolap_noinv = wfolap_noinv + &
290 291 292
                           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
293 294
                  lm = lm + nn
               END DO
295
            END DO
Matthias Redies's avatar
Matthias Redies committed
296
         END DO
297
      END DO
Matthias Redies's avatar
Matthias Redies committed
298 299

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

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

!       wfolap = wfolap + cdum

Matthias Redies's avatar
Matthias Redies committed
310 311
   END FUNCTION wfolap_noinv
END MODULE m_olap