abcof.F90 16.1 KB
Newer Older
1 2
MODULE m_abcof
CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
3
  SUBROUTINE abcof(input,atoms,sym, cell,lapw,ne,usdus,&
4 5
       noco,jspin,oneD, acof,bcof,ccof,zMat,&
       eig,acoflo,bcoflo,e1cof,e2cof,aveccof,bveccof,cveccof)
6 7 8 9 10 11 12
    !     ************************************************************
    !     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
Daniel Wortmann's avatar
Daniel Wortmann committed
13
    USE m_setabc1lo
14 15 16 17 18 19
    USE m_sphbes
    USE m_dsphbs
    USE m_abclocdn
    USE m_ylm
    USE m_types
    IMPLICIT NONE
20 21
    TYPE(t_input),INTENT(IN)  :: input
    TYPE(t_usdus),INTENT(IN)  :: usdus
22 23 24 25 26 27
    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
28
    TYPE(t_zMat),INTENT(IN)   :: zMat
29 30
    !     ..
    !     .. Scalar Arguments ..
Markus Betzinger's avatar
Markus Betzinger committed
31 32
    INTEGER, INTENT (IN) :: ne
    INTEGER, INTENT (IN) :: jspin
33 34
    !     ..
    !     .. Array Arguments ..
35 36 37 38 39 40 41 42 43 44 45
    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)
    REAL,    OPTIONAL, INTENT (IN)  :: eig(:)!(dimension%neigd)
    COMPLEX, OPTIONAL, INTENT (OUT) :: acoflo(-atoms%llod:,:,:,:)
    COMPLEX, OPTIONAL, INTENT (OUT) :: bcoflo(-atoms%llod:,:,:,:)
    COMPLEX, OPTIONAL, INTENT (OUT) :: e1cof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, OPTIONAL, INTENT (OUT) :: e2cof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, OPTIONAL, INTENT (OUT) :: aveccof(:,:,0:,:)!(3,nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, OPTIONAL, INTENT (OUT) :: bveccof(:,:,0:,:)!(3,nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, OPTIONAL, INTENT (OUT) :: cveccof(:,-atoms%llod:,:,:,:)
46 47 48
    !     ..
    !     .. Local Scalars ..
    COMPLEX cexp,phase,c_0,c_1,c_2,ci
Daniel Wortmann's avatar
Daniel Wortmann committed
49
    REAL const,df,r1,s,tmk,wronk,qss(3)
50
    REAL s2h, s2h_e(ne)
Daniel Wortmann's avatar
Daniel Wortmann committed
51
    INTEGER i,j,k,l,ll1,lm ,n,nap,natom,nn,iatom,jatom,lmp,m,nkvec
52 53 54
    INTEGER inv_f,ie,ilo,kspin,iintsp,nintsp,nvmax,lo,inap
    !     ..
    !     .. Local Arrays ..
Daniel Wortmann's avatar
Daniel Wortmann committed
55
    INTEGER nbasf0(atoms%nlod,atoms%nat)
56
    REAL dfj(0:atoms%lmaxd),fj(0:atoms%lmaxd),fg(3),fgp(3),fgr(3),fk(3),fkp(3),fkr(3)
Daniel Wortmann's avatar
Daniel Wortmann committed
57
    REAL alo1(atoms%nlod),blo1(atoms%nlod),clo1(atoms%nlod)
58 59
    COMPLEX ylm( (atoms%lmaxd+1)**2 )
    COMPLEX ccchi(2,2)
Daniel Wortmann's avatar
Daniel Wortmann committed
60
    LOGICAL enough(atoms%nat),apw(0:atoms%lmaxd,atoms%ntype)
61 62
    REAL,    ALLOCATABLE :: work_r(:)
    COMPLEX, ALLOCATABLE :: work_c(:)
63

64 65
    
    IF (zmat%l_real) THEN
66 67 68 69
       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

70
    !     ..
71 72
    ci = CMPLX(0.0,1.0)
    const = 2 * tpi_const/SQRT(cell%omtil)
73
    !
74 75 76 77 78 79 80 81 82 83 84 85
    acof(:,:,:)   = CMPLX(0.0,0.0)
    bcof(:,:,:)   = CMPLX(0.0,0.0)
    ccof(:,:,:,:) = CMPLX(0.0,0.0)
    IF(PRESENT(eig)) THEN
       acoflo(:,:,:,:)    = CMPLX(0.0,0.0)
       bcoflo(:,:,:,:)    = CMPLX(0.0,0.0)
       e1cof(:,:,:)       = CMPLX(0.0,0.0)
       e2cof(:,:,:)       = CMPLX(0.0,0.0)
       aveccof(:,:,:,:)   = CMPLX(0.0,0.0)
       bveccof(:,:,:,:)   = CMPLX(0.0,0.0)
       cveccof(:,:,:,:,:) = CMPLX(0.0,0.0)
    END IF
86 87 88 89
    !     ..
    !+APW_LO
    DO n = 1, atoms%ntype
       DO l = 0,atoms%lmax(n)
90
          apw(l,n) = .FALSE.
91
          DO lo = 1,atoms%nlo(n)
92
             IF (atoms%l_dulo(lo,n)) apw(l,n) = .TRUE.
93
          ENDDO
94
          IF ((input%l_useapw).AND.(atoms%lapw_l(n).GE.l)) apw(l,n) = .FALSE.
95 96
       ENDDO
       DO lo = 1,atoms%nlo(n)
97
          IF (atoms%l_dulo(lo,n)) apw(atoms%llo(lo,n),n) = .TRUE.
98 99 100 101 102 103 104 105 106 107 108 109
       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)
       IF (iintsp .EQ. 1) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
110
          qss= - noco%qss/2
111
       ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
112
          qss= + noco%qss/2
113 114 115
       ENDIF

       !---> loop over atom types
116 117
       !$OMP PARALLEL DO &
       !$OMP& DEFAULT(none)&
118 119 120
       !$OMP& PRIVATE(n,nn,natom,k,i,work_r,work_c,ccchi,kspin,fg,fk,s,r1,fj,dfj,l,df,wronk,tmk,phase,&
       !$OMP& alo1,blo1,clo1,inap,nap,j,fgr,fgp,s2h,s2h_e,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,eig,&
Daniel Wortmann's avatar
Daniel Wortmann committed
121 122 123
       !$OMP& jspin,qss,&
       !$OMP& apw,const,&
       !$OMP& nbasf0,enough,&
124
       !$OMP& acof,bcof,ccof,e1cof,e2cof,acoflo,bcoflo,aveccof,bveccof,cveccof)
125
       DO n = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
126 127
          CALL setabc1lo(atoms,n,usdus,jspin,alo1,blo1,clo1)
          
128 129
          !  ----> loop over equivalent atoms
          DO nn = 1,atoms%neq(n)
130 131 132 133 134
             natom = 0
             DO i = 1, n-1
                natom = natom + atoms%neq(i)
             ENDDO
             natom = natom + nn
135 136
             IF ((atoms%invsat(natom).EQ.0) .OR. (atoms%invsat(natom).EQ.1)) THEN
                !--->        loop over lapws
137
                IF (zmat%l_real) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
138
                   ALLOCATE ( work_r(ne) )
139
                ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
140
                   ALLOCATE ( work_c(ne) )
141
                ENDIF
142 143
                DO k = 1,nvmax
                   IF (.NOT.noco%l_noco) THEN
144
                      IF (zmat%l_real) THEN
145
                         work_r(:ne)=zMat%z_r(k,:ne)
146
                      ELSE
147
                         work_c(:ne)=zMat%z_c(k,:ne)
148
                      END IF
149 150 151 152
                   ENDIF

                   IF (noco%l_noco) THEN
                      !--->            generate the complex conjgates of the spinors (chi)
153 154 155 156
                      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))
157 158 159 160 161
                      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
162
                            work_c(i) = ccchi(iintsp,jspin)*zMat%z_c(kspin+k,i)
163 164 165 166 167 168 169
                         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
170 171
                            work_c(i) = ccchi(1,jspin)*zMat%z_c(k,i)&
                                 &                        + ccchi(2,jspin)*zMat%z_c(kspin+k,i)
172 173 174 175
                         ENDDO
                      ENDIF
                   ENDIF ! (noco%l_noco)
                   IF (noco%l_ss) THEN
176
                      fg = lapw%gvec(:,k,iintsp) + qss
177
                   ELSE
178
                      fg = lapw%gvec(:,k,jspin) + qss
179
                   ENDIF ! (noco%l_ss)
180 181 182 183 184 185
                   fk = lapw%bkpt + fg
                   s =  DOT_PRODUCT(fk,MATMUL(cell%bbmat,fk))
                   IF(PRESENT(eig)) THEN
                      s2h = 0.5*s
                      s2h_e(:ne) = s2h-eig(:ne)
                   END IF
186
                   s = SQRT(s)
187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
                   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)+&
203 204
                                     fk(2)*atoms%taual(2,natom)+&
                                     fk(3)*atoms%taual(3,natom))
205
                   phase = CMPLX(COS(tmk),SIN(tmk))
206 207 208 209 210 211 212
                   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
213 214
                      fkr(j) = 0.0
                      fgr(j) = 0.0
215 216 217
                      DO i = 1,3
                         IF (oneD%odi%d1) THEN
                            fkr(j) = fkr(j) + fk(i)*oneD%ods%mrot(i,j,inap)
218
                            fgr(j) = fgr(j) + fg(i)*oneD%ods%mrot(i,j,inap)
219 220
                         ELSE
                            fkr(j) = fkr(j) + fk(i)*sym%mrot(i,j,inap)
221
                            fgr(j) = fgr(j) + fg(i)*sym%mrot(i,j,inap)
222
                         END IF
223 224 225 226
                      END DO
                   END DO
                   fkp = MATMUL(fkr,cell%bmat)
                   fgp = MATMUL(fgr,cell%bmat)
227 228 229 230 231 232 233 234
                   !     ----> 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
235
                         c_0 = CONJG(ylm(lm+1))*phase
236 237 238
                         c_1 = c_0 *  fj(l)
                         c_2 = c_0 * dfj(l)
                         !     ----> loop over bands
239
                         IF (zmat%l_real) THEN
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
                            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)
256 257 258 259 260 261 262 263
                               IF ((atoms%l_geo(n)).AND.(PRESENT(eig))) THEN
                                  CALL CPP_BLAS_caxpy(ne,c_1,work_c*s2h_e,1, e1cof(1,lmp,jatom),1)
                                  CALL CPP_BLAS_caxpy(ne,c_2,work_c*s2h_e,1, e2cof(1,lmp,jatom),1)
                                  DO i = 1,3
                                     CALL CPP_BLAS_caxpy(ne,c_1,work_c*fgp(i),1, aveccof(i,1,lmp,jatom),3)
                                     CALL CPP_BLAS_caxpy(ne,c_2,work_c*fgp(i),1, bveccof(i,1,lmp,jatom),3)
                                  END DO
                               END IF
264
                            ENDIF
265 266 267
                         ENDIF
                      ENDDO ! loop over m
                   ENDDO ! loop over l
Daniel Wortmann's avatar
Daniel Wortmann committed
268 269 270 271
                   DO lo=1,atoms%nlo(n)
                      DO nkvec=1,lapw%nkvec(lo,natom)
                         IF (k==lapw%kvec(nkvec,lo,natom)) THEN !check if this k-vector has LO attached
                            CALL abclocdn(atoms,sym,noco,lapw,cell,ccchi(:,jspin),iintsp,phase,ylm,&
272 273
                                 n,natom,k,nkvec,lo,ne,alo1,blo1,clo1,acof,bcof,ccof,zMat,&
                                 fgp,acoflo,bcoflo,aveccof,bveccof,cveccof)
Daniel Wortmann's avatar
Daniel Wortmann committed
274 275 276
                         ENDIF
                      ENDDO
                   END DO
277
                ENDDO ! loop over LAPWs
278
                IF (zmat%l_real) THEN
279 280 281 282
                   DEALLOCATE(work_r)
                ELSE               
                   DEALLOCATE(work_c)
                ENDIF
283 284 285
             ENDIF  ! invsatom == ( 0 v 1 )
          ENDDO    ! loop over equivalent atoms
       ENDDO       ! loop over atom types
286
       !$OMP END PARALLEL DO
287
    ENDDO       ! loop over interstitial spin
288
    IF (noco%l_soc.AND.sym%invs) THEN
289

290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
       !
       !                           -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)&
Daniel Wortmann's avatar
Daniel Wortmann committed
311
                     &             + atoms%taual(:,iatom),lapw%bkpt))
312 313 314 315 316
                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
317 318 319 320 321 322 323 324
                         ccof(m,ie,ilo,jatom) = inv_f * cexp * CONJG(  ccof(-m,ie,ilo,iatom))
                         IF(PRESENT(eig)) THEN
                            acoflo(m,ie,ilo,jatom) = inv_f * cexp * CONJG(acoflo(-m,ie,ilo,iatom))
                            bcoflo(m,ie,ilo,jatom) = inv_f * cexp * CONJG(bcoflo(-m,ie,ilo,iatom))
                            DO i = 1,3
                               cveccof(i,m,ie,ilo,jatom) = -inv_f * cexp * CONJG(cveccof(i,-m,ie,ilo,iatom))
                            END DO
                         END IF
325
                      ENDDO
326 327
                   ENDDO
                ENDDO
328 329 330 331 332 333 334 335 336
                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))
                         bcof(ie,lm,jatom) = inv_f * cexp * CONJG(bcof(ie,lmp,iatom))
337 338 339 340 341 342 343 344 345 346 347
                      END DO
                      IF ((atoms%l_geo(n)).AND.(PRESENT(eig))) THEN
                         DO ie = 1,ne
                            e1cof(ie,lm,jatom) = inv_f * cexp * CONJG(e1cof(ie,lmp,iatom))
                            e2cof(ie,lm,jatom) = inv_f * cexp * CONJG(e2cof(ie,lmp,iatom))
                            DO i = 1,3
                               aveccof(i,ie,lm,jatom) = -inv_f * cexp * CONJG(aveccof(i,ie,lmp,iatom))
                               bveccof(i,ie,lm,jatom) = -inv_f * cexp * CONJG(bveccof(i,ie,lmp,iatom))
                            END DO
                         END DO
                      END IF
348 349
                   ENDDO
                ENDDO
350 351
             ENDIF
          ENDDO
352
       ENDDO
353
    ENDIF
354 355 356

  END SUBROUTINE abcof
END MODULE m_abcof