abcof.F90 13.5 KB
Newer Older
1 2
MODULE m_abcof
CONTAINS
3
  SUBROUTINE abcof(input,atoms,nobd,sym, cell, bkpt,lapw,ne,usdus,&
4
       noco,jspin,kveclo,oneD, acof,bcof,ccof,zMat)
5 6 7 8 9 10 11 12 13 14 15 16 17 18
    !     ************************************************************
    !     subroutine constructs the a,b coefficients of the linearized
    !     m.t. wavefunctions for each band and atom.       c.l. fu
    !     ************************************************************
#include "cpp_double.h"

    USE m_constants, ONLY : tpi_const
    USE m_setabc1locdn
    USE m_sphbes
    USE m_dsphbs
    USE m_abclocdn
    USE m_ylm
    USE m_types
    IMPLICIT NONE
19 20
    TYPE(t_input),INTENT(IN)  :: input
    TYPE(t_usdus),INTENT(IN)  :: usdus
21 22 23 24 25 26
    TYPE(t_lapw),INTENT(IN)   :: lapw
    TYPE(t_oneD),INTENT(IN)   :: oneD
    TYPE(t_noco),INTENT(IN)   :: noco
    TYPE(t_sym),INTENT(IN)    :: sym
    TYPE(t_cell),INTENT(IN)   :: cell
    TYPE(t_atoms),INTENT(IN)  :: atoms
27
    TYPE(t_zMat),INTENT(IN)   :: zMat
28 29
    !     ..
    !     .. Scalar Arguments ..
Markus Betzinger's avatar
Markus Betzinger committed
30 31 32
    INTEGER, INTENT (IN) :: nobd
    INTEGER, INTENT (IN) :: ne
    INTEGER, INTENT (IN) :: jspin
33 34 35 36
    !     ..
    !     .. Array Arguments ..
    INTEGER, INTENT (IN) :: kveclo(atoms%nlotot)
    REAL,    INTENT (IN) :: bkpt(3)
Daniel Wortmann's avatar
Daniel Wortmann committed
37 38 39
    COMPLEX, INTENT (OUT):: acof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, INTENT (OUT):: bcof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, INTENT (OUT):: ccof(-atoms%llod:,:,:,:)!(-llod:llod,nobd,atoms%nlod,atoms%nat)
40 41 42 43 44 45 46 47
    !     ..
    !     .. Local Scalars ..
    COMPLEX cexp,phase,c_0,c_1,c_2,ci
    REAL const,df,r1,s,tmk,wronk,qss1,qss2,qss3
    INTEGER i,j,k,l,ll1,lm ,n,nap,natom,nn,iatom,jatom,lmp,m
    INTEGER inv_f,ie,ilo,kspin,iintsp,nintsp,nvmax,lo,inap
    !     ..
    !     .. Local Arrays ..
Daniel Wortmann's avatar
Daniel Wortmann committed
48 49
    INTEGER kvec(2*(2*atoms%llod+1),atoms%nlod,atoms%nat  )
    INTEGER nbasf0(atoms%nlod,atoms%nat),nkvec(atoms%nlod,atoms%nat)
50
    REAL dfj(0:atoms%lmaxd),fj(0:atoms%lmaxd),fk(3),fkp(3),fkr(3)
Daniel Wortmann's avatar
Daniel Wortmann committed
51
    REAL alo1(atoms%nlod,atoms%ntype),blo1(atoms%nlod,atoms%ntype),clo1(atoms%nlod,atoms%ntype)
52 53
    COMPLEX ylm( (atoms%lmaxd+1)**2 )
    COMPLEX ccchi(2,2)
Daniel Wortmann's avatar
Daniel Wortmann committed
54
    LOGICAL enough(atoms%nat),apw(0:atoms%lmaxd,atoms%ntype)
55 56
    REAL,    ALLOCATABLE :: work_r(:)
    COMPLEX, ALLOCATABLE :: work_c(:)
57

58 59
    
    IF (zmat%l_real) THEN
60 61 62 63
       IF (noco%l_soc.AND.sym%invs) CALL judft_error("BUG in abcof, SOC&INVS but real?")
       IF (noco%l_noco) CALL judft_error("BUG in abcof, l_noco but real?")
    ENDIF

64
    !     ..
65 66
    ci = CMPLX(0.0,1.0)
    const = 2 * tpi_const/SQRT(cell%omtil)
67
    !
68 69
    acof(:,:,:) = CMPLX(0.0,0.0)
    bcof(:,:,:) = CMPLX(0.0,0.0)
70 71 72 73
    !     ..
    !+APW_LO
    DO n = 1, atoms%ntype
       DO l = 0,atoms%lmax(n)
74
          apw(l,n) = .FALSE.
75
          DO lo = 1,atoms%nlo(n)
76
             IF (atoms%l_dulo(lo,n)) apw(l,n) = .TRUE.
77
          ENDDO
78
          IF ((input%l_useapw).AND.(atoms%lapw_l(n).GE.l)) apw(l,n) = .FALSE.
79 80
       ENDDO
       DO lo = 1,atoms%nlo(n)
81
          IF (atoms%l_dulo(lo,n)) apw(atoms%llo(lo,n),n) = .TRUE.
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
       ENDDO
    ENDDO
    !+APW_LO
    !
    nintsp = 1
    IF (noco%l_ss) nintsp = 2
    !---> loop over the interstitial spin
    DO iintsp = 1,nintsp
       !
       nvmax=lapw%nv(jspin)
       IF (noco%l_ss) nvmax=lapw%nv(iintsp)
       CALL setabc1locdn(jspin, atoms,lapw,ne,noco,iintsp, sym,usdus,&
            kveclo, enough,nkvec,kvec,nbasf0,ccof, alo1,blo1,clo1)
       !
       IF (iintsp .EQ. 1) THEN
          qss1= - noco%qss(1)/2
          qss2= - noco%qss(2)/2
          qss3= - noco%qss(3)/2
       ELSE
          qss1= + noco%qss(1)/2
          qss2= + noco%qss(2)/2
          qss3= + noco%qss(3)/2
       ENDIF

       !---> loop over atom types
107 108 109 110 111 112 113 114 115
       !$OMP PARALLEL DO &
       !$OMP& DEFAULT(none)&
       !$OMP& PRIVATE(n,nn,natom,k,i,work_r,work_c,ccchi,kspin,fk,s,r1,fj,dfj,l,df,wronk,tmk,phase,&
       !$OMP& inap,nap,j,fkr,fkp,ylm,ll1,m,c_0,c_1,c_2,jatom,lmp,inv_f,lm)&
       !$OMP& SHARED(noco,atoms,sym,cell,oneD,lapw,nvmax,ne,zMat,usdus,ci,iintsp,&
       !$OMP& jspin,bkpt,qss1,qss2,qss3,&
       !$OMP& apw,const,nobd,&
       !$OMP& alo1,blo1,clo1,kvec,nbasf0,nkvec,enough,&
       !$OMP& acof,bcof,ccof)
116 117 118
       DO n = 1,atoms%ntype
          !  ----> loop over equivalent atoms
          DO nn = 1,atoms%neq(n)
119 120 121 122 123
             natom = 0
             DO i = 1, n-1
                natom = natom + atoms%neq(i)
             ENDDO
             natom = natom + nn
124 125
             IF ((atoms%invsat(natom).EQ.0) .OR. (atoms%invsat(natom).EQ.1)) THEN
                !--->        loop over lapws
126
                IF (zmat%l_real) THEN
127 128 129 130
                   ALLOCATE ( work_r(nobd) )
                ELSE
                   ALLOCATE ( work_c(nobd) )
                ENDIF
131 132
                DO k = 1,nvmax
                   IF (.NOT.noco%l_noco) THEN
133
                      IF (zmat%l_real) THEN
134
                         work_r(:ne)=zMat%z_r(k,:ne)
135
                      ELSE
136
                         work_c(:ne)=zMat%z_c(k,:ne)
137
                      END IF
138 139 140 141
                   ENDIF

                   IF (noco%l_noco) THEN
                      !--->            generate the complex conjgates of the spinors (chi)
142 143 144 145
                      ccchi(1,1) = CONJG( EXP(-ci*noco%alph(n)/2)*COS(noco%beta(n)/2))
                      ccchi(1,2) = CONJG(-EXP(-ci*noco%alph(n)/2)*SIN(noco%beta(n)/2))
                      ccchi(2,1) = CONJG( EXP( ci*noco%alph(n)/2)*SIN(noco%beta(n)/2))
                      ccchi(2,2) = CONJG( EXP( ci*noco%alph(n)/2)*COS(noco%beta(n)/2))
146 147 148 149 150
                      IF (noco%l_ss) THEN
                         !--->              the coefficients of the spin-down basis functions are
                         !--->              stored in the second half of the eigenvector
                         kspin = (iintsp-1)*(lapw%nv(1)+atoms%nlotot)
                         DO i = 1,ne
151
                            work_c(i) = ccchi(iintsp,jspin)*zMat%z_c(kspin+k,i)
152 153 154 155 156 157 158
                         ENDDO
                      ELSE
                         !--->              perform sum over the two interstitial spin directions
                         !--->              and take into account the spin boundary conditions
                         !--->              (jspin counts the local spin directions inside each MT)
                         kspin = lapw%nv(1)+atoms%nlotot
                         DO i = 1,ne
159 160
                            work_c(i) = ccchi(1,jspin)*zMat%z_c(k,i)&
                                 &                        + ccchi(2,jspin)*zMat%z_c(kspin+k,i)
161 162 163 164 165 166 167 168 169 170 171 172
                         ENDDO
                      ENDIF
                   ENDIF ! (noco%l_noco)
                   IF (noco%l_ss) THEN
                      fk(1) = bkpt(1) + lapw%k1(k,iintsp) + qss1
                      fk(2) = bkpt(2) + lapw%k2(k,iintsp) + qss2
                      fk(3) = bkpt(3) + lapw%k3(k,iintsp) + qss3
                   ELSE
                      fk(1) = bkpt(1) + lapw%k1(k,jspin) + qss1
                      fk(2) = bkpt(2) + lapw%k2(k,jspin) + qss2
                      fk(3) = bkpt(3) + lapw%k3(k,jspin) + qss3
                   ENDIF ! (noco%l_ss)
173 174
                   s=  DOT_PRODUCT(fk,MATMUL(cell%bbmat,fk))
                   s = SQRT(s)
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192
                   r1 = atoms%rmt(n)*s
                   CALL sphbes(atoms%lmax(n),r1,fj)
                   CALL dsphbs(atoms%lmax(n),r1,fj,dfj)
                   !  ----> construct a and b coefficients
                   DO l = 0,atoms%lmax(n)
                      df = s*dfj(l)
                      wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin) - usdus%us(l,n,jspin)*usdus%duds(l,n,jspin)
                      IF (apw(l,n)) THEN
                         fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin)
                         dfj(l) = 0.0d0
                      ELSE
                         dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk
                         fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk
                      ENDIF
                   ENDDO ! loop over l
                   tmk = tpi_const* (fk(1)*atoms%taual(1,natom)+&
                        &                     fk(2)*atoms%taual(2,natom)+&
                        &                     fk(3)*atoms%taual(3,natom))
193
                   phase = CMPLX(COS(tmk),SIN(tmk))
194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209
                   IF (oneD%odi%d1) THEN
                      inap = oneD%ods%ngopr(natom)
                   ELSE
                      nap = atoms%ngopr(natom)
                      inap = sym%invtab(nap)
                   END IF
                   DO j = 1,3
                      fkr(j) = 0.
                      DO i = 1,3
                         IF (oneD%odi%d1) THEN
                            fkr(j) = fkr(j) + fk(i)*oneD%ods%mrot(i,j,inap)
                         ELSE
                            fkr(j) = fkr(j) + fk(i)*sym%mrot(i,j,inap)
                         END IF
                      ENDDO
                   ENDDO
210
                   fkp=MATMUL(fkr,cell%bmat)
211 212 213 214 215 216 217 218
                   !     ----> generate spherical harmonics
                   CALL ylm4(atoms%lmax(n),fkp,ylm)
                   !     ----> loop over l
                   DO l = 0,atoms%lmax(n)
                      ll1 = l* (l+1)
                      !     ----> loop over m
                      DO m = -l,l
                         lm = ll1 + m
219
                         c_0 = CONJG(ylm(lm+1))*phase
220 221 222
                         c_1 = c_0 *  fj(l)
                         c_2 = c_0 * dfj(l)
                         !     ----> loop over bands
223
                         IF (zmat%l_real) THEN
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
                            acof(:ne,lm,natom) = acof(:ne,lm,natom) +  c_1 * work_r(:ne)
                            bcof(:ne,lm,natom) = bcof(:ne,lm,natom) +  c_2 * work_r(:ne)
                         ELSE
                            acof(:ne,lm,natom) = acof(:ne,lm,natom) +  c_1 * work_c(:ne)
                            bcof(:ne,lm,natom) = bcof(:ne,lm,natom) +  c_2 * work_c(:ne)
                         END IF

                         IF (noco%l_soc.AND.sym%invs) THEN
                            IF (atoms%invsat(natom).EQ.1) THEN
                               jatom = sym%invsatnr(natom)
                               lmp = ll1 - m
                               inv_f = (-1)**(l-m)
                               c_1 =  CONJG(c_1) * inv_f
                               c_2 =  CONJG(c_2) * inv_f
                               CALL CPP_BLAS_caxpy(ne,c_1,work_c,1, acof(1,lmp,jatom),1)
                               CALL CPP_BLAS_caxpy(ne,c_2,work_c,1, bcof(1,lmp,jatom),1)
                            ENDIF
241 242 243 244 245
                         ENDIF
                      ENDDO ! loop over m
                   ENDDO ! loop over l
                   IF (.NOT.enough(natom)) THEN
                      CALL abclocdn(atoms,sym, noco,ccchi(1,jspin),kspin,iintsp,const,phase,ylm,n,natom,k,&
246
                           s,nvmax,ne,nbasf0,alo1,blo1,clo1,kvec(1,1,natom),nkvec,enough(natom),acof,bcof,ccof,zMat)
247 248
                   ENDIF
                ENDDO ! loop over LAPWs
249
                IF (zmat%l_real) THEN
250 251 252 253
                   DEALLOCATE(work_r)
                ELSE               
                   DEALLOCATE(work_c)
                ENDIF
254 255 256
             ENDIF  ! invsatom == ( 0 v 1 )
          ENDDO    ! loop over equivalent atoms
       ENDDO       ! loop over atom types
257
       !$OMP END PARALLEL DO
258
    ENDDO       ! loop over interstitial spin
259
    IF (noco%l_soc.AND.sym%invs) THEN
260

261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
       !
       !                           -p,n       (l+m)   p,n  *
       ! Usually, we exploit that A     = (-1)      (A    )  if p and -p are the positions
       !                           l,m                l,-m
       ! of two atoms related by inversion symmetry and the coefficients are considered to
       ! be in the local frame of the representative atom. This is possible, if z is real.
       ! After SOC, however, the eigenvectors z are complex and this is no longer possible
       ! so the z has to enter, not z*. This is done within the k-loop.
       !                                    -p,n       m   p,n  *
       ! When called from hsohelp, we need A     = (-1)  (A    ) because we don't have to
       !                                     l,m           l,-m                    l
       ! rotate, but in the sums in hsoham only products A*  A   enter and the (-1) cancels.
       !                                                  lm  lm
    ELSE
       iatom = 0
       DO n = 1,atoms%ntype
          DO nn = 1,atoms%neq(n)
             iatom = iatom + 1
             IF (atoms%invsat(iatom).EQ.1) THEN
                jatom = sym%invsatnr(iatom)
                cexp = EXP(tpi_const*ci*DOT_PRODUCT(atoms%taual(:,jatom)&
                     &             + atoms%taual(:,iatom),(/bkpt(1),bkpt(2),bkpt(3)/)))
                DO ilo = 1,atoms%nlo(n)
                   l = atoms%llo(ilo,n)
                   DO m = -l,l
                      inv_f = (-1.0)**(m+l)
                      DO ie = 1,ne
                         ccof(m,ie,ilo,jatom) = inv_f * cexp *CONJG(  ccof(-m,ie,ilo,iatom))
                      ENDDO
290 291
                   ENDDO
                ENDDO
292 293 294 295 296 297 298 299 300 301 302 303
                DO l = 0,atoms%lmax(n)
                   ll1 = l* (l+1)
                   DO m =-l,l
                      lm  = ll1 + m
                      lmp = ll1 - m
                      inv_f = (-1.0)**(m+l)
                      DO ie = 1,ne
                         acof(ie,lm,jatom) = inv_f * cexp * CONJG(acof(ie,lmp,iatom))
                      ENDDO
                      DO ie = 1,ne
                         bcof(ie,lm,jatom) = inv_f * cexp * CONJG(bcof(ie,lmp,iatom))
                      ENDDO
304 305
                   ENDDO
                ENDDO
306 307
             ENDIF
          ENDDO
308
       ENDDO
309
    ENDIF
310 311 312

  END SUBROUTINE abcof
END MODULE m_abcof