hlomat.F90 16.6 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13
MODULE m_hlomat
!***********************************************************************
! updates the hamiltonian  matrix with the contributions from the local
! orbitals.
! p.kurz sept. 1996
!***********************************************************************
CONTAINS
14 15
  SUBROUTINE hlomat(input,atoms,mpi,lapw,ud,tlmplm,sym,cell,noco,isp,&
       ntyp,na,fj,gj,alo1,blo1,clo1, iintsp,jintsp,chi,hmat)
16
    !
17
    USE m_hsmt_ab
18 19
    USE m_types
    IMPLICIT NONE
20
    TYPE(t_input),INTENT(IN)  :: input
21 22
    TYPE(t_atoms),INTENT(IN)  :: atoms
    TYPE(t_lapw),INTENT(IN)   :: lapw
23 24 25 26 27 28
    TYPE(t_mpi),INTENT(IN)    :: mpi
    TYPE(t_usdus),INTENT(IN)  :: ud
    TYPE(t_tlmplm),INTENT(IN) :: tlmplm
    TYPE(t_sym),INTENT(IN)    :: sym
    TYPE(t_cell),INTENT(IN)   :: cell
    TYPE(t_noco),INTENT(IN)   :: noco
29 30


31 32 33
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: na,ntyp  
34
    INTEGER, INTENT (IN) :: isp !spin for usdus and tlmplm
35
    INTEGER, INTENT (IN) :: jintsp,iintsp
36
    COMPLEX, INTENT (IN) :: chi
37 38
    !     ..
    !     .. Array Arguments ..
39 40
    REAL, INTENT (IN) :: alo1(:),blo1(:),clo1(:)
    REAL,INTENT(IN)      :: fj(:,0:,:),gj(:,0:,:)
41

42
    CLASS(t_mat),INTENT (INOUT) :: hmat
43 44
    !     ..
    !     .. Local Scalars ..
45 46 47 48
    COMPLEX axx,bxx,cxx,dtd,dtu,dtulo,ulotd,ulotu,ulotulo,utd,utu, utulo
    INTEGER im,in,invsfct,l,lm,lmp,lo,lolo,lolop,lop,lp,i  
    INTEGER mp,nkvec,nkvecp,lmplm,loplo,kp,m,mlo,mlolo
    INTEGER locol,lorow,ii,ij,n,k,ab_size
49 50
    !     ..
    !     .. Local Arrays ..
51
    COMPLEX, ALLOCATABLE :: ab(:,:,:),ax(:),bx(:),cx(:)
52
    COMPLEX,ALLOCATABLE  :: abclo(:,:,:,:,:)
53
    !     ..
54

55

56
    !-->              synthesize the complex conjugates of a and b
57
    ALLOCATE(ab(MAXVAL(lapw%nv),0:2*atoms%lmaxd*(atoms%lmaxd+2)+1,MIN(jintsp,iintsp):MAX(jintsp,iintsp)))
58 59
    ALLOCATE(ax(MAXVAL(lapw%nv)),bx(MAXVAL(lapw%nv)),cx(MAXVAL(lapw%nv)))
    ALLOCATE(abclo(3,-atoms%llod:atoms%llod,2*(2*atoms%llod+1),atoms%nlod,2))
60
    DO i=MIN(jintsp,iintsp),MAX(jintsp,iintsp)
61
       CALL hsmt_ab(sym,atoms,noco,isp,i,ntyp,na,cell,lapw,fj,gj,ab(:,:,i),ab_size,.TRUE.,abclo(:,:,:,:,i),alo1,blo1,clo1)
62
    ENDDO
63

64

65 66 67 68 69 70 71 72
    mlo=0;mlolo=0
    DO m=1,ntyp-1
       mlo=mlo+atoms%nlo(m)
       mlolo=mlolo+atoms%nlo(m)*(atoms%nlo(m)+1)/2
    ENDDO


    !$OMP MASTER
Matthias Redies's avatar
Matthias Redies committed
73
    IF ((atoms%invsat(na) == 0) .OR. (atoms%invsat(na) == 1)) THEN
74 75 76 77 78 79
       !--->    if this atom is the first of two atoms related by inversion,
       !--->    the contributions to the overlap matrix of both atoms are added
       !--->    at once. where it is made use of the fact, that the sum of
       !--->    these contributions is twice the real part of the contribution
       !--->    of each atom. note, that in this case there are twice as many
       !--->    (2*(2*l+1)) k-vectors (compare abccoflo and comments there).
Matthias Redies's avatar
Matthias Redies committed
80 81
       IF (atoms%invsat(na) == 0) invsfct = 1
       IF (atoms%invsat(na) == 1) invsfct = 2
82
       !
83

84 85 86 87 88 89
       DO lo = 1,atoms%nlo(ntyp)
          l = atoms%llo(lo,ntyp)
          !--->       calculate the hamiltonian matrix elements with the regular
          !--->       flapw basis-functions
          DO m = -l,l
             lm = l* (l+1) + m
90
             DO kp = 1,lapw%nv(iintsp)
91 92 93
                ax(kp) = CMPLX(0.0,0.0)
                bx(kp) = CMPLX(0.0,0.0)
                cx(kp) = CMPLX(0.0,0.0)
94
             END DO
95
             DO lp = 0,atoms%lnonsph(ntyp)
96 97
                DO mp = -lp,lp
                   lmp = lp* (lp+1) + mp
98
                   in = tlmplm%ind(lmp,lm,ntyp,isp)
99
                   IF (lmp==lm) in=(lm* (lm+3))/2
100 101
                   IF (in.NE.-9999) THEN
                      IF (in.GE.0) THEN
102 103 104 105
                         utu = tlmplm%tuu(in,ntyp,isp)
                         dtu = tlmplm%tdu(in,ntyp,isp)
                         utd = tlmplm%tud(in,ntyp,isp)
                         dtd = tlmplm%tdd(in,ntyp,isp)
106 107
                      ELSE
                         im = -in
108 109 110 111
                         utu = CONJG(tlmplm%tuu(im,ntyp,isp))
                         dtu = CONJG(tlmplm%tud(im,ntyp,isp))
                         utd = CONJG(tlmplm%tdu(im,ntyp,isp))
                         dtd = CONJG(tlmplm%tdd(im,ntyp,isp))
112
                      END IF
113 114
                      utulo = tlmplm%tuulo(lmp,m,lo+mlo,isp)
                      dtulo = tlmplm%tdulo(lmp,m,lo+mlo,isp)
115 116 117 118
                      !--->                   note, that utu,dtu... are the t-matrices and
                      !--->                   not their complex conjugates as in hssphn
                      !--->                   and that a,b,alo... are the complex
                      !--->                   conjugates of the a,b...-coefficients
119
                      DO kp = 1,lapw%nv(iintsp)
120 121 122
                         ax(kp) = ax(kp) + ab(kp,lmp,iintsp)*utu + ab(kp,ab_size/2+lmp,iintsp)*dtu
                         bx(kp) = bx(kp) + ab(kp,lmp,iintsp)*utd + ab(kp,ab_size/2+lmp,iintsp)*dtd
                         cx(kp) = cx(kp) + ab(kp,lmp,iintsp)*utulo + ab(kp,ab_size/2+lmp,iintsp)*dtulo
123 124 125 126 127 128
                      END DO
                   END IF
                END DO
             END DO
             !+t3e
             DO nkvec = 1,invsfct* (2*l+1)
129 130
                locol= lapw%nv(jintsp)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix
                IF (MOD(locol-1,mpi%n_size) == mpi%n_rank) THEN !only this MPI rank calculates this column
131 132
                   locol=(locol-1)/mpi%n_size+1 !this is the column in local storage
                   IF (hmat%l_real) THEN
133
                      DO kp = 1,lapw%nv(iintsp)
134
                         hmat%data_r(kp,locol) = hmat%data_r(kp,locol) + chi*invsfct * (&
135 136 137 138 139 140
                              REAL(abclo(1,m,nkvec,lo,jintsp))* REAL(ax(kp)) +&
                              AIMAG(abclo(1,m,nkvec,lo,jintsp))*AIMAG(ax(kp)) +&
                              REAL(abclo(2,m,nkvec,lo,jintsp))* REAL(bx(kp)) +&
                              AIMAG(abclo(2,m,nkvec,lo,jintsp))*AIMAG(bx(kp)) +&
                              REAL(abclo(3,m,nkvec,lo,jintsp))* REAL(cx(kp)) +&
                              AIMAG(abclo(3,m,nkvec,lo,jintsp))*AIMAG(cx(kp)) )
141 142
                         IF (input%l_useapw) THEN
                            !---> APWlo
143
                            hmat%data_r(kp,locol) = hmat%data_r(kp,locol) + 0.25 * atoms%rmt(ntyp)**2 * chi*invsfct * (&
144 145
                                 (CONJG(ab(kp,lm,iintsp))* ud%us(l,ntyp,isp)+&
                                 CONJG(ab(kp,ab_size/2+lm,iintsp))*ud%uds(l,ntyp,isp))*&
146 147 148
                                 (abclo(1,m,nkvec,lo,jintsp)*  ud%dus(l,ntyp,isp)&
                                 +abclo(2,m,nkvec,lo,jintsp)* ud%duds(l,ntyp,isp)&
                                 +abclo(3,m,nkvec,lo,jintsp)*ud%dulos(lo,ntyp,isp) ))
149
                         ENDIF
150 151
                      ENDDO
                   ELSE
152
                      DO kp = 1,lapw%nv(iintsp)
153 154 155 156 157 158 159 160 161 162 163 164 165
                         hmat%data_c(kp,locol) = hmat%data_c(kp,locol) + chi*invsfct * (&
                              abclo(1,m,nkvec,lo,jintsp) * CONJG( ax(kp) ) +&
                              abclo(2,m,nkvec,lo,jintsp) * CONJG( bx(kp) ) +&
                              abclo(3,m,nkvec,lo,jintsp) * CONJG( cx(kp) ) )
                         IF (input%l_useapw) THEN
                            !---> APWlo
                            hmat%data_c(kp,locol)=hmat%data_c(kp,locol) + 0.25 * atoms%rmt(ntyp)**2 * chi*invsfct*(&
                                 (CONJG(ab(kp,lm,iintsp))* ud%us(l,ntyp,isp)+&
                                 CONJG(ab(kp,ab_size/2+lm,iintsp))*ud%uds(l,ntyp,isp))*&
                                 (abclo(1,m,nkvec,lo,jintsp)*  ud%dus(l,ntyp,isp)&
                                 +abclo(2,m,nkvec,lo,jintsp)* ud%duds(l,ntyp,isp)&
                                 +abclo(3,m,nkvec,lo,jintsp)*ud%dulos(lo,ntyp,isp) ))
                         ENDIF
166 167
                      ENDDO
                   ENDIF
168 169 170 171 172 173 174
                   !--->             jump to the last matrixelement of the current row
                ENDIF
             END DO
          END DO
          !--->       calculate the hamiltonian matrix elements with other
          !--->       local orbitals at the same atom and with itself
          DO nkvec = 1,invsfct* (2*l+1)
175 176
             locol= lapw%nv(jintsp)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix
             IF (MOD(locol-1,mpi%n_size) == mpi%n_rank) THEN !only this MPI rank calculates this column
177
                locol=(locol-1)/mpi%n_size+1 !this is the column in local storage
178 179
                !--->          calculate the hamiltonian matrix elements with other
                !--->          local orbitals at the same atom, if they have the same l
180 181
                DO lop = 1, MERGE(lo-1,atoms%nlo(ntyp),iintsp==jintsp)
                   IF (lop==lo) CYCLE
182 183
                   lp = atoms%llo(lop,ntyp)
                   DO nkvecp = 1,invsfct* (2*lp+1)
184
                      lorow=lapw%nv(iintsp)+lapw%index_lo(lop,na)+nkvecp
185 186 187 188
                      DO m = -l,l
                         lm = l* (l+1) + m
                         DO mp = -lp,lp
                            lmp = lp* (lp+1) + mp
189
                            in = tlmplm%ind(lmp,lm,ntyp,isp)
190
                            IF (lmp==lm) in=(lm* (lm+3))/2
191 192
                            IF (in.NE.-9999) THEN
                               IF (in.GE.0) THEN
193 194 195 196
                                  utu = tlmplm%tuu(in,ntyp,isp)
                                  dtu = tlmplm%tdu(in,ntyp,isp)
                                  utd = tlmplm%tud(in,ntyp,isp)
                                  dtd = tlmplm%tdd(in,ntyp,isp)
197 198
                               ELSE
                                  im = -in
199 200 201 202
                                  utu = CONJG(tlmplm%tuu(im,ntyp,isp))
                                  dtu = CONJG(tlmplm%tud(im,ntyp,isp))
                                  utd = CONJG(tlmplm%tdu(im,ntyp,isp))
                                  dtd = CONJG(tlmplm%tdd(im,ntyp,isp))
203
                               END IF
204 205 206 207
                               utulo = tlmplm%tuulo(lmp,m,lo+mlo,isp)
                               dtulo = tlmplm%tdulo(lmp,m,lo+mlo,isp)
                               ulotu=CONJG(tlmplm%tuulo(lm,mp,lop+mlo,isp))
                               ulotd=CONJG(tlmplm%tdulo(lm,mp,lop+mlo,isp))
208
                               !--->                         note that lo > lop
209 210 211 212 213 214 215
                               IF (lo>lop) THEN
                                  lolop = ((lo-1)*lo)/2 + lop
                                  ulotulo = CONJG(tlmplm%tuloulo (m,mp,lolop+mlolo,isp))
                               ELSE
                                  lolop = ((lop-1)*lop)/2 + lo
                                  ulotulo = CONJG(tlmplm%tuloulo (mp,m,lolop+mlolo,isp))
                               ENDIF
216 217 218 219 220 221
                               axx=CONJG(abclo(1,m,nkvec,lo,jintsp))*utu +&
                                    CONJG(abclo(2,m,nkvec,lo,jintsp))*utd +&
                                    CONJG(abclo(3,m,nkvec,lo,jintsp))*utulo
                               bxx=CONJG(abclo(1,m,nkvec,lo,jintsp))*dtu +&
                                    CONJG(abclo(2,m,nkvec,lo,jintsp))*dtd +&
                                    CONJG(abclo(3,m,nkvec,lo,jintsp))*dtulo
222
                               cxx = &
223 224 225
                                    CONJG(abclo(1,m,nkvec,lo,jintsp))*ulotu +&
                                    CONJG(abclo(2,m,nkvec,lo,jintsp))*ulotd +&
                                    CONJG(abclo(3,m,nkvec,lo,jintsp))*ulotulo
226 227
                               IF (hmat%l_real) THEN
                                  hmat%data_r(lorow,locol) = hmat%data_r(lorow,locol) + chi*invsfct * (&
228 229 230 231 232 233
                                       REAL(abclo(1,mp,nkvecp,lop,iintsp))* REAL(axx) -&
                                       AIMAG(abclo(1,mp,nkvecp,lop,iintsp))*AIMAG(axx) +&
                                       REAL(abclo(2,mp,nkvecp,lop,iintsp))* REAL(bxx) -&
                                       AIMAG(abclo(2,mp,nkvecp,lop,iintsp))*AIMAG(bxx) +&
                                       REAL(abclo(3,mp,nkvecp,lop,iintsp))* REAL(cxx) -&
                                       AIMAG(abclo(3,mp,nkvecp,lop,iintsp))*AIMAG(cxx) )
234
                               ELSE
235
                                  hmat%data_c(lorow,locol) = hmat%data_c(lorow,locol) + chi*invsfct * CONJG(&
236 237 238
                                       abclo(1,mp,nkvecp,lop,iintsp) * axx +&
                                       abclo(2,mp,nkvecp,lop,iintsp) * bxx +&
                                       abclo(3,mp,nkvecp,lop,iintsp) * cxx )
239 240 241 242 243 244 245 246
                               ENDIF
                            END IF
                         END DO
                      END DO
                   END DO
                END DO
                !--->          calculate the hamiltonian matrix elements of one local
                !--->          orbital with itself
247 248
                lop=lo
                DO nkvecp = 1,MERGE(nkvec,invsfct* (2*l+1),iintsp==jintsp)
249
                   lorow=lapw%nv(iintsp)+lapw%index_lo(lop,na)+nkvecp
250 251 252 253
                   DO m = -l,l
                      lm = l* (l+1) + m
                      DO mp = -l,l
                         lmp = l* (l+1) + mp
254
                         in = tlmplm%ind(lmp,lm,ntyp,isp)
255
                         IF (lmp==lm) in=(lm* (lm+3))/2
256 257
                         IF (in.NE.-9999) THEN
                            IF (in.GE.0) THEN
258 259 260 261
                               utu = tlmplm%tuu(in,ntyp,isp)
                               dtu = tlmplm%tdu(in,ntyp,isp)
                               utd = tlmplm%tud(in,ntyp,isp)
                               dtd = tlmplm%tdd(in,ntyp,isp)
262 263
                            ELSE
                               im = -in
264 265 266 267
                               utu = CONJG(tlmplm%tuu(im,ntyp,isp))
                               dtu = CONJG(tlmplm%tud(im,ntyp,isp))
                               utd = CONJG(tlmplm%tdu(im,ntyp,isp))
                               dtd = CONJG(tlmplm%tdd(im,ntyp,isp))
268
                            END IF
269 270 271 272
                            utulo = tlmplm%tuulo(lmp,m,lo+mlo,isp)
                            dtulo = tlmplm%tdulo(lmp,m,lo+mlo,isp)
                            ulotu = CONJG(tlmplm%tuulo(lm,mp,lo+mlo,isp))
                            ulotd = CONJG(tlmplm%tdulo(lm,mp,lo+mlo,isp))
273
                            lolo = ((lo-1)*lo)/2 + lo
274
                            ulotulo =CONJG(tlmplm%tuloulo(m,mp,lolo+mlolo,isp))
275 276 277 278 279 280 281 282 283
                            axx = CONJG(abclo(1,m,nkvec,lo,jintsp))*utu +&
                                 CONJG(abclo(2,m,nkvec,lo,jintsp))*utd +&
                                 CONJG(abclo(3,m,nkvec,lo,jintsp))*utulo
                            bxx = CONJG(abclo(1,m,nkvec,lo,jintsp))*dtu +&
                                 CONJG(abclo(2,m,nkvec,lo,jintsp))*dtd +&
                                 CONJG(abclo(3,m,nkvec,lo,jintsp))*dtulo
                            cxx = CONJG(abclo(1,m,nkvec,lo,jintsp))*ulotu +&
                                 CONJG(abclo(2,m,nkvec,lo,jintsp))*ulotd +&
                                 CONJG(abclo(3,m,nkvec,lo,jintsp))*ulotulo
284 285
                            IF (hmat%l_real) THEN
                               hmat%data_r(lorow,locol) = hmat%data_r(lorow,locol) + chi*invsfct* (&
286 287 288 289 290 291
                                    REAL(abclo(1,mp,nkvecp,lo,iintsp))* REAL(axx) -&
                                    AIMAG(abclo(1,mp,nkvecp,lo,iintsp))*AIMAG(axx) +&
                                    REAL(abclo(2,mp,nkvecp,lo,iintsp))* REAL(bxx) -&
                                    AIMAG(abclo(2,mp,nkvecp,lo,iintsp))*AIMAG(bxx) +&
                                    REAL(abclo(3,mp,nkvecp,lo,iintsp))* REAL(cxx) -&
                                    AIMAG(abclo(3,mp,nkvecp,lo,iintsp))*AIMAG(cxx) )
292
                            ELSE
293
                               hmat%data_c(lorow,locol) = hmat%data_c(lorow,locol) + chi*invsfct* CONJG(&
294 295 296
                                    abclo(1,mp,nkvecp,lo,iintsp)*axx +&
                                    abclo(2,mp,nkvecp,lo,iintsp)*bxx +&
                                    abclo(3,mp,nkvecp,lo,iintsp)*cxx )
297 298 299 300 301
                            ENDIF
                         END IF
                      END DO
                   END DO
                END DO
302
             ENDIF !If this lo to be calculated by mpi rank
303 304
          END DO
       END DO ! end of lo = 1,atoms%nlo loop
305

306 307 308 309
    END IF
    !$OMP END MASTER
    !$OMP barrier
  END SUBROUTINE hlomat
310
END MODULE m_hlomat