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

11
    USE m_constants, ONLY : tpi_const, ImagUnit
12
    USE m_setabc1lo
13 14 15 16 17
    USE m_sphbes
    USE m_dsphbs
    USE m_abclocdn
    USE m_ylm
    USE m_types
18
    USE m_juDFT
19
    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_mat),INTENT(IN)    :: zMat
29
    TYPE(t_force),OPTIONAL,INTENT(INOUT) :: force
30 31
    !     ..
    !     .. Scalar Arguments ..
Markus Betzinger's avatar
Markus Betzinger committed
32 33
    INTEGER, INTENT (IN) :: ne
    INTEGER, INTENT (IN) :: jspin
34 35
    !     ..
    !     .. Array Arguments ..
36 37 38
    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)
39
    REAL,    OPTIONAL, INTENT (IN) :: eig(:)!(dimension%neigd)
40 41
    !     ..
    !     .. Local Scalars ..
42
    COMPLEX cexp,phase,c_0,c_1,c_2
43
    REAL const,df,r1,s,tmk,wronk,qss(3)
44
    REAL s2h, s2h_e(ne)
45
    INTEGER i,j,k,l,ll1,lm ,n,nap,natom,nn,iatom,jatom,lmp,m,nkvec
46
    INTEGER inv_f,ie,ilo,kspin,iintsp,nintsp,nvmax,lo,inap
47
    LOGICAL l_force
48 49
    !     ..
    !     .. Local Arrays ..
50
    INTEGER nbasf0(atoms%nlod,atoms%nat)
51
    REAL dfj(0:atoms%lmaxd),fj(0:atoms%lmaxd),fg(3),fgp(3),fgr(3),fk(3),fkp(3),fkr(3)
52
    REAL alo1(atoms%nlod),blo1(atoms%nlod),clo1(atoms%nlod)
53 54
    COMPLEX ylm( (atoms%lmaxd+1)**2 )
    COMPLEX ccchi(2,2)
55
    LOGICAL enough(atoms%nat),apw(0:atoms%lmaxd,atoms%ntype)
56 57
    REAL,    ALLOCATABLE :: work_r(:)
    COMPLEX, ALLOCATABLE :: work_c(:)
58

59 60
    CALL timestart("abcof")

61
    IF (zmat%l_real) THEN
62 63 64 65 66
       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

    const = 2 * tpi_const/SQRT(cell%omtil)
67

68 69 70
    acof(:,:,:)   = CMPLX(0.0,0.0)
    bcof(:,:,:)   = CMPLX(0.0,0.0)
    ccof(:,:,:,:) = CMPLX(0.0,0.0)
71 72 73 74 75
    l_force = .FALSE.
    IF(PRESENT(eig).AND.input%l_f) THEN
       l_force = .TRUE.
    END IF
    IF(l_force) THEN
76 77 78 79 80 81 82
       force%acoflo  = CMPLX(0.0,0.0)
       force%bcoflo  = CMPLX(0.0,0.0)
       force%e1cof   = CMPLX(0.0,0.0)
       force%e2cof   = CMPLX(0.0,0.0)
       force%aveccof = CMPLX(0.0,0.0)
       force%bveccof = CMPLX(0.0,0.0)
       force%cveccof = CMPLX(0.0,0.0)
83
    END IF
84

85 86 87
    !+APW_LO
    DO n = 1, atoms%ntype
       DO l = 0,atoms%lmax(n)
88
          apw(l,n) = .FALSE.
89
          DO lo = 1,atoms%nlo(n)
90
             IF (atoms%l_dulo(lo,n)) apw(l,n) = .TRUE.
91
          ENDDO
92
          IF ((input%l_useapw).AND.(atoms%lapw_l(n).GE.l)) apw(l,n) = .FALSE.
93 94
       ENDDO
       DO lo = 1,atoms%nlo(n)
95
          IF (atoms%l_dulo(lo,n)) apw(atoms%llo(lo,n),n) = .TRUE.
96 97 98 99 100 101 102 103 104 105 106 107
       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
108
          qss= - noco%qss/2
109
       ELSE
110
          qss= + noco%qss/2
111 112 113
       ENDIF

       !---> loop over atom types
114 115
       !$OMP PARALLEL DO &
       !$OMP& DEFAULT(none)&
116
       !$OMP& PRIVATE(n,nn,natom,k,i,work_r,work_c,ccchi,kspin,fg,fk,s,r1,fj,dfj,l,df,wronk,tmk,phase,lo,nkvec,&
117
       !$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)&
118
       !$OMP& SHARED(noco,atoms,sym,cell,oneD,lapw,nvmax,ne,zMat,usdus,iintsp,eig,l_force,&
119
       !$OMP& jspin,qss,apw,const,nbasf0,enough,acof,bcof,ccof,force)
120
       DO n = 1,atoms%ntype
121 122
          CALL setabc1lo(atoms,n,usdus,jspin,alo1,blo1,clo1)
          
123 124
          !  ----> loop over equivalent atoms
          DO nn = 1,atoms%neq(n)
125 126 127 128 129
             natom = 0
             DO i = 1, n-1
                natom = natom + atoms%neq(i)
             ENDDO
             natom = natom + nn
130 131
             IF ((atoms%invsat(natom).EQ.0) .OR. (atoms%invsat(natom).EQ.1)) THEN
                !--->        loop over lapws
132
                IF (zmat%l_real) THEN
133
                   ALLOCATE ( work_r(ne) )
134
                ELSE
135
                   ALLOCATE ( work_c(ne) )
136
                ENDIF
137 138
                DO k = 1,nvmax
                   IF (.NOT.noco%l_noco) THEN
139
                      IF (zmat%l_real) THEN
140
                         work_r(:ne)=zMat%data_r(k,:ne)
141
                      ELSE
142
                         work_c(:ne)=zMat%data_c(k,:ne)
143
                      END IF
144 145 146
                   ENDIF

                   IF (noco%l_noco) THEN
147 148 149 150 151
                      !--->            generate the spinors (chi)
                      ccchi(1,1) =  EXP(ImagUnit*noco%alph(n)/2)*COS(noco%beta(n)/2)
                      ccchi(1,2) = -EXP(ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2)
                      ccchi(2,1) =  EXP(-ImagUnit*noco%alph(n)/2)*SIN(noco%beta(n)/2)
                      ccchi(2,2) =  EXP(-ImagUnit*noco%alph(n)/2)*COS(noco%beta(n)/2)
152 153 154 155 156
                      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
157
                            work_c(i) = ccchi(iintsp,jspin)*zMat%data_c(kspin+k,i)
158 159 160 161 162 163 164
                         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
165 166
                            work_c(i) = ccchi(1,jspin)*zMat%data_c(k,i)&
                                 &                        + ccchi(2,jspin)*zMat%data_c(kspin+k,i)
167 168 169 170
                         ENDDO
                      ENDIF
                   ENDIF ! (noco%l_noco)
                   IF (noco%l_ss) THEN
171
                      fg = lapw%gvec(:,k,iintsp) + qss
172
                   ELSE
173
                      fg = lapw%gvec(:,k,jspin) + qss
174
                   ENDIF ! (noco%l_ss)
175 176
                   fk = lapw%bkpt + fg
                   s =  DOT_PRODUCT(fk,MATMUL(cell%bbmat,fk))
177
                   IF(l_force) THEN
178 179 180
                      s2h = 0.5*s
                      s2h_e(:ne) = s2h-eig(:ne)
                   END IF
181
                   s = SQRT(s)
182 183 184 185 186 187 188 189 190
                   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)
Matthias Redies's avatar
Matthias Redies committed
191
                         dfj(l) = 0.0
192 193 194 195 196 197
                      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)+&
198 199
                                     fk(2)*atoms%taual(2,natom)+&
                                     fk(3)*atoms%taual(3,natom))
200
                   phase = CMPLX(COS(tmk),SIN(tmk))
201 202 203 204 205 206 207
                   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
208 209
                      fkr(j) = 0.0
                      fgr(j) = 0.0
210 211 212
                      DO i = 1,3
                         IF (oneD%odi%d1) THEN
                            fkr(j) = fkr(j) + fk(i)*oneD%ods%mrot(i,j,inap)
213
                            fgr(j) = fgr(j) + fg(i)*oneD%ods%mrot(i,j,inap)
214 215
                         ELSE
                            fkr(j) = fkr(j) + fk(i)*sym%mrot(i,j,inap)
216
                            fgr(j) = fgr(j) + fg(i)*sym%mrot(i,j,inap)
217
                         END IF
218 219 220 221
                      END DO
                   END DO
                   fkp = MATMUL(fkr,cell%bmat)
                   fgp = MATMUL(fgr,cell%bmat)
222 223 224 225 226 227 228 229
                   !     ----> 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
230
                         c_0 = CONJG(ylm(lm+1))*phase
231 232 233
                         c_1 = c_0 *  fj(l)
                         c_2 = c_0 * dfj(l)
                         !     ----> loop over bands
234
                         IF (zmat%l_real) THEN
235 236 237 238 239 240 241
                            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

242
                         IF (atoms%l_geo(n).AND.l_force) THEN
243
                            IF (zmat%l_real) THEN
244 245
                               force%e1cof(:ne,lm,natom) = force%e1cof(:ne,lm,natom) + c_1 * work_r(:ne) * s2h_e(:ne)
                               force%e2cof(:ne,lm,natom) = force%e2cof(:ne,lm,natom) + c_2 * work_r(:ne) * s2h_e(:ne)
246
                               DO i = 1,3
247 248
                                  force%aveccof(i,:ne,lm,natom) = force%aveccof(i,:ne,lm,natom) + c_1 * work_r(:ne) * fgp(i)
                                  force%bveccof(i,:ne,lm,natom) = force%bveccof(i,:ne,lm,natom) + c_2 * work_r(:ne) * fgp(i)
249 250
                               END DO
                            ELSE
251 252
                               force%e1cof(:ne,lm,natom) = force%e1cof(:ne,lm,natom) + c_1 * work_c(:ne) * s2h_e(:ne)
                               force%e2cof(:ne,lm,natom) = force%e2cof(:ne,lm,natom) + c_2 * work_c(:ne) * s2h_e(:ne)
253
                               DO i = 1,3
254 255
                                  force%aveccof(i,:ne,lm,natom) = force%aveccof(i,:ne,lm,natom) + c_1 * work_c(:ne) * fgp(i)
                                  force%bveccof(i,:ne,lm,natom) = force%bveccof(i,:ne,lm,natom) + c_2 * work_c(:ne) * fgp(i)
256 257 258 259
                               END DO
                            END IF
                         ENDIF

260 261 262 263 264 265 266 267 268
                         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)
269
                               IF (atoms%l_geo(n).AND.l_force) THEN
270 271
                                  CALL CPP_BLAS_caxpy(ne,c_1,work_c*s2h_e,1, force%e1cof(1,lmp,jatom),1)
                                  CALL CPP_BLAS_caxpy(ne,c_2,work_c*s2h_e,1, force%e2cof(1,lmp,jatom),1)
272
                                  DO i = 1,3
273 274
                                     CALL CPP_BLAS_caxpy(ne,c_1,work_c*fgp(i),1, force%aveccof(i,1,lmp,jatom),3)
                                     CALL CPP_BLAS_caxpy(ne,c_2,work_c*fgp(i),1, force%bveccof(i,1,lmp,jatom),3)
275 276
                                  END DO
                               END IF
277
                            ENDIF
278 279 280
                         ENDIF
                      ENDDO ! loop over m
                   ENDDO ! loop over l
281 282 283 284
                   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,&
285
                                          n,natom,k,nkvec,lo,ne,alo1,blo1,clo1,acof,bcof,ccof,zMat,l_force,fgp,force)
286 287 288
                         ENDIF
                      ENDDO
                   END DO
289
                ENDDO ! loop over LAPWs
290
                IF (zmat%l_real) THEN
291 292 293 294
                   DEALLOCATE(work_r)
                ELSE               
                   DEALLOCATE(work_c)
                ENDIF
295 296 297
             ENDIF  ! invsatom == ( 0 v 1 )
          ENDDO    ! loop over equivalent atoms
       ENDDO       ! loop over atom types
298
       !$OMP END PARALLEL DO
299
    ENDDO       ! loop over interstitial spin
300
    IF (noco%l_soc.AND.sym%invs) THEN
301

302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
       !
       !                           -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)
322
                cexp = EXP(tpi_const*ImagUnit*DOT_PRODUCT(atoms%taual(:,jatom)&
323
                     &             + atoms%taual(:,iatom),lapw%bkpt))
324 325 326 327 328
                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
329
                         ccof(m,ie,ilo,jatom) = inv_f * cexp * CONJG(  ccof(-m,ie,ilo,iatom))
330
                         IF(l_force) THEN
331 332
                            force%acoflo(m,ie,ilo,jatom) = inv_f * cexp * CONJG(force%acoflo(-m,ie,ilo,iatom))
                            force%bcoflo(m,ie,ilo,jatom) = inv_f * cexp * CONJG(force%bcoflo(-m,ie,ilo,iatom))
333
                            DO i = 1,3
334
                               force%cveccof(i,m,ie,ilo,jatom) = -inv_f * cexp * CONJG(force%cveccof(i,-m,ie,ilo,iatom))
335 336
                            END DO
                         END IF
337
                      ENDDO
338 339
                   ENDDO
                ENDDO
340 341 342 343 344 345 346 347 348
                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))
349
                      END DO
350
                      IF (atoms%l_geo(n).AND.l_force) THEN
351
                         DO ie = 1,ne
352 353
                            force%e1cof(ie,lm,jatom) = inv_f * cexp * CONJG(force%e1cof(ie,lmp,iatom))
                            force%e2cof(ie,lm,jatom) = inv_f * cexp * CONJG(force%e2cof(ie,lmp,iatom))
354
                            DO i = 1,3
355 356
                               force%aveccof(i,ie,lm,jatom) = -inv_f * cexp * CONJG(force%aveccof(i,ie,lmp,iatom))
                               force%bveccof(i,ie,lm,jatom) = -inv_f * cexp * CONJG(force%bveccof(i,ie,lmp,iatom))
357 358 359
                            END DO
                         END DO
                      END IF
360 361
                   ENDDO
                ENDDO
362 363
             ENDIF
          ENDDO
364
       ENDDO
365
    ENDIF
366

367 368
    CALL timestop("abcof")

369 370
  END SUBROUTINE abcof
END MODULE m_abcof