locrectify.F90 19.8 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 14 15
!*************************************************************
!     Find linear combinations of local orbitals that are
!     eigenfunctions of the z-reflection operation.
!     Used to exploit z-reflection symmetry when solving
!     the eigenvalue problem (film calculations).
!     Frank Freimuth, January 2007
!*************************************************************
#include"cpp_double.h"
MODULE m_locrectify
16
  USE m_juDFT
17
  CONTAINS
18 19
    SUBROUTINE locrectify(jsp,input,lapw,bkpt,atoms, kveclo, sym,cell,&
         kindlocrec,evenlocs,oddlocs,evenlocindex, oddlocindex,locrec_r,locrec_c)
20

21 22 23 24 25 26 27 28 29 30 31 32 33
      USE m_ylm
      USE m_constants,ONLY:tpi_const
      USE m_loccoeff
      USE m_types
      IMPLICIT NONE
      INTEGER,INTENT(IN)         :: jsp
      TYPE(t_input),INTENT(IN)   :: input
      TYPE(t_sym),INTENT(IN)     :: sym
      TYPE(t_cell),INTENT(IN)    :: cell
      TYPE(t_atoms),INTENT(IN)   :: atoms
      TYPE(t_lapw),INTENT(IN)    :: lapw
      REAL,INTENT(in)::bkpt(3)
      INTEGER,INTENT(in)::kveclo(atoms%nlotot)
34

35 36
      REAL,OPTIONAL,INTENT(out)::locrec_r(atoms%nlotot,atoms%nlotot)
      COMPLEX,OPTIONAL,INTENT(out)::locrec_c(atoms%nlotot,atoms%nlotot)
37

38 39 40 41 42
      LOGICAL,INTENT(out)::kindlocrec(atoms%nlotot)
      INTEGER,INTENT(out)::evenlocs
      INTEGER,INTENT(out)::oddlocs
      INTEGER,INTENT(out)::evenlocindex(atoms%nlotot)
      INTEGER,INTENT(out)::oddlocindex(atoms%nlotot)
43

44 45 46 47 48 49 50 51 52 53 54 55 56
      INTEGER numorec
      INTEGER invatom,zrefatom
      LOGICAL l_planetwo,l_spacetwo,l_spacefour,l_spacetwoinv
      LOGICAL l_crosstwo
      INTEGER pos,loc,nn,angmom,bas,locm,loopi,j,i
      INTEGER lm,locmm,lmm
      INTEGER  ind,vec,natom,nap,inap,locnum,ll1
      INTEGER ipiv(28),info
      COMPLEX loccoffs(28,28),reccoffs(28)
      COMPLEX ylm(16),phase
      REAL fk(3),tmk,fkp(3)
      !      complex coeffi(1:28,1:28,0:13)     transferred to mod_loccoeff.f
      !      logical coffkind(1:28,0:13)        import with USE m_loccoeff
Daniel Wortmann's avatar
Daniel Wortmann committed
57 58
      INTEGER jump(atoms%nat)
      INTEGER mqnum,nalo,naup,basindex(atoms%nat,atoms%nlod)
59 60
      INTEGER zrefnap,zrefinap,coffindex
      LOGICAL l_real,l_invsatswap
61

62
      !      ci=cmplx(0.0,1.0) transferred to mod_loccoeff.f
63

64
      l_real=PRESENT(locrec_r)
65

Daniel Wortmann's avatar
Daniel Wortmann committed
66
      jump(1:atoms%nat)=0
67 68 69 70 71 72
      IF (l_real) THEN
         locrec_r=0.0
      ELSE
         locrec_c=0.0
      ENDIF
      numorec=0
73

74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
      !***********************************************
      !   basindex holds first index of subset of locs
      !***********************************************
      bas=1
      nalo=1
      DO ind=1,atoms%ntype
         naup=nalo+atoms%neq(ind)-1
         DO natom=nalo,naup
            IF(atoms%invsat(natom).EQ.2)CYCLE
            DO loc=1,atoms%nlo(ind)
               basindex(natom,loc)=bas
               locnum=(atoms%invsat(natom)+1)*(2*atoms%llo(loc,ind)+1)
               bas=bas+locnum
            ENDDO !loc
         ENDDO !natom
         nalo=naup+1
      ENDDO !ind
91

92 93 94 95 96 97 98 99 100 101
      bas=0
      nalo=1
      DO ind=1,atoms%ntype
         naup=nalo+atoms%neq(ind)-1
         DO natom=nalo,naup
            IF(atoms%invsat(natom).EQ.2)CYCLE
            IF(jump(natom).NE.0)THEN
               bas=bas+jump(natom)
               CYCLE
            ENDIF
102

103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
            !*******************************************************
            !       find out by what kind of symmetry the atoms
            !       are interrelated
            !*******************************************************
            !two atoms related by inversion symmetry are in xy-plane
            l_planetwo=.FALSE.    
            !two atoms are related by z-reflection, but not by inversion
            l_spacetwo=.FALSE.
            !two atoms are related by z-reflection and by inversion
            l_spacetwoinv=.FALSE.
            !four atoms are entangled, by both inversion and z-reflection
            l_spacefour=.FALSE.
            !two atoms are related by inversion, but not by z-reflection
            l_crosstwo=.FALSE.
            !order of locs corresponds to setup of coeffi-array => l_invsatswap=.false.
            !order of locs is reversed with respect to coeffi-array => l_invsatswap=.true.
            l_invsatswap=.FALSE.
            IF((ABS(ABS(atoms%taual(3,natom))-0.5).LT.1e-6).AND..NOT.input%film)THEN !atom on face
               IF(atoms%invsat(natom).EQ.1)THEN ! two atoms
                  invatom=sym%invsatnr(natom)
                  l_planetwo=.TRUE.
               ENDIF
            ELSEIF(ABS(atoms%taual(3,natom)).GT.1e-6) THEN !atom not in xy-plane
               DO zrefatom=nalo,naup !find z-reflection image
                  IF(ALL(ABS(atoms%taual(1:2,natom)-atoms%taual(1:2,zrefatom))-&
                       NINT(ABS(atoms%taual(1:2,natom)-atoms%taual(1:2,zrefatom))).LT.1e-6)) THEN
                     IF(ABS(atoms%taual(3,natom)+atoms%taual(3,zrefatom)).LT.1e-6) GOTO 543
                  ENDIF
               ENDDO !zrefatom
               if (l_real) THEN
                  !if there is no z-reflected counterpart there must be an inversion 
                  !counterpart instead
                  IF(atoms%invsat(natom).EQ.1)THEN
                     zrefatom=0 !this switches on l_crosstwo later on
                  ELSE
                     CALL juDFT_error("nozref-invsat",calledby ="locrectify")
                  ENDIF
               else
                  !in the complex version there must be a z-reflected counterpart
                  CALL juDFT_error("zrefatom not found",calledby ="locrectify")
               endif
543            CONTINUE
145

146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161
               IF(atoms%invsat(natom).EQ.0) THEN !no inversion-image
                  l_spacetwo=.TRUE.
                  if (l_real)   CALL juDFT_error("spacetwo & INVERSION",calledby ="locrectify")
               ELSE                        !inversion-image
                  invatom=sym%invsatnr(natom)
                  IF(zrefatom.EQ.0)THEN
                     l_crosstwo=.TRUE.
                  ELSEIF(invatom.EQ.zrefatom)THEN    !inversion = z-reflection
                     l_spacetwoinv=.TRUE.
                  ELSE                           !inversion /= z-reflection
                     if (l_real) THEN
                        l_spacefour=.TRUE. !four entangled locs
                     else
                        CALL juDFT_error("spacefour & no INVERSION",calledby ="locrectify")
                     endif
                     IF(atoms%invsat(zrefatom).EQ.1)THEN
162

163 164
                     ELSEIF(atoms%invsat(zrefatom).EQ.2)THEN
                        l_invsatswap=.TRUE.
Daniel Wortmann's avatar
Daniel Wortmann committed
165
                        DO info=1,atoms%nat
166 167 168 169 170 171 172 173 174 175
                           IF(sym%invsatnr(info).EQ.zrefatom)GOTO 723
                        ENDDO
                        CALL juDFT_error("no atoms for zref",calledby ="locrectify")
723                     CONTINUE
                        zrefatom=info
                        IF(atoms%invsat(zrefatom)/=1) CALL juDFT_error("atoms%invsat" ,calledby ="locrectify")
                        IF(.NOT.ALL(ABS(atoms%taual(1:2,invatom)- atoms%taual(1:2,zrefatom))<1e-6)) &
                             CALL juDFT_error ("invsat-zref2",calledby ="locrectify")
                        IF(.NOT.ABS(atoms%taual(3,invatom)+atoms%taual(3,zrefatom)) <1e-6)  &
                             CALL juDFT_error("invsat-zref3", calledby="locrectify")
176

177 178 179 180 181 182 183 184 185 186 187
                     ELSEIF(atoms%invsat(zrefatom).EQ.0)THEN
                        CALL juDFT_error("spacefour: atoms-zref",calledby ="locrectify")
                     ENDIF
                  ENDIF
               ENDIF
            ELSE !atom lies in the xy-plane
               IF(atoms%invsat(natom).EQ.1)THEN ! two atoms
                  invatom=sym%invsatnr(natom)
                  l_planetwo=.TRUE.
               ENDIF
            ENDIF !symmetry partners
188

189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282
            DO loc=1,atoms%nlo(ind)
               angmom=atoms%llo(loc,ind)
               IF(angmom.EQ.0)THEN
                  IF(l_planetwo)THEN
                     kindlocrec(numorec+1)=.TRUE.
                     kindlocrec(numorec+2)=.TRUE.
                     IF (l_real) THEN
                        locrec_r(bas+1,numorec+1)=1.0
                        locrec_r(bas+2,numorec+2)=1.0
                     ELSE
                        locrec_c(bas+1,numorec+1)=1.0
                        locrec_c(bas+2,numorec+2)=1.0
                     ENDIF
                     bas=bas+2
                     numorec=numorec+2
                     CYCLE
                  ELSEIF(l_spacetwoinv.OR.l_crosstwo)THEN
                     coffindex=0
                     locnum=2
                  ELSEIF(l_spacetwo) THEN
                     locnum=1
                     coffindex=0
                  ELSEIF(l_spacefour)THEN
                     locnum=2
                     coffindex=6
                  ELSE
                     kindlocrec(numorec+1)=.TRUE.
                     IF (l_real) THEN
                        locrec_r(bas+1,numorec+1)=1.0
                     ELSE
                        locrec_c(bas+1,numorec+1)=1.0
                     END IF
                     bas=bas+1
                     numorec=numorec+1
                     CYCLE
                  ENDIF !l_planetwo,l_crosstwo,....
               ELSE
                  IF(l_planetwo)THEN
                     locnum=2*(2*angmom+1)
                     IF(angmom.EQ.1)THEN
                        coffindex=3
                     ELSEIF(angmom.EQ.2)THEN
                        coffindex=9
                     ELSEIF(angmom.EQ.3)THEN
                        coffindex=12
                     ELSE
                        CALL juDFT_error("angmom1",calledby ="locrectify")
                     ENDIF
                  ELSEIF(l_spacetwoinv.OR.l_crosstwo)THEN
                     locnum=2*(2*angmom+1)
                     IF(angmom.EQ.1)THEN
                        coffindex=4
                     ELSEIF(angmom.EQ.2)THEN
                        coffindex=8
                     ELSEIF(angmom.EQ.3)THEN
                        coffindex=11
                     ELSE
                        CALL juDFT_error("angmom5",calledby ="locrectify")
                     ENDIF
                  ELSEIF(l_spacetwo) THEN
                     locnum=2*angmom+1
                     IF(angmom.EQ.1)THEN
                        coffindex=4
                     ELSEIF(angmom.EQ.2)THEN
                        coffindex=8
                     ELSEIF(angmom.EQ.3)THEN
                        coffindex=11
                     ELSE
                        CALL juDFT_error("angmom2",calledby="locrectify")
                     ENDIF
                  ELSEIF(l_spacefour) THEN
                     locnum=2*(2*angmom+1)
                     IF(angmom.EQ.1)THEN
                        coffindex=7
                     ELSEIF(angmom.EQ.2)THEN
                        coffindex=10
                     ELSEIF(angmom.EQ.3)THEN
                        coffindex=13
                     ELSE
                        CALL juDFT_error("angmom3",calledby ="locrectify")
                     ENDIF
                  ELSE
                     locnum=2*angmom+1
                     IF(angmom.EQ.1)THEN
                        coffindex=1
                     ELSEIF(angmom.EQ.2)THEN
                        coffindex=2
                     ELSEIF(angmom.EQ.3)THEN
                        coffindex=5
                     ELSE
                        CALL juDFT_error("angmom",calledby ="locrectify")
                     ENDIF
                  ENDIF! l_planetwo
               ENDIF !angmom
283

284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
               mqnum=2*angmom+1
               ll1=angmom*angmom-1
               loccoffs(:,:)=CMPLX(0.0,0.0)
               !**************************************************
               !        Write the expansion of the locs
               !        in terms of spherical harmonics into the
               !        loccoffs-array.
               !        First index of loccoffs: m-quantum-number.
               !        Second index of loccoffs: index of loc.
               !**************************************************
               DO locm=1,locnum
                  pos=bas+locm
                  vec=kveclo(pos)
                  !The vector of the planewave onto which the loc is matched.
                  fk(:)=bkpt(:)+(/lapw%k1(vec,jsp),lapw%k2(vec,jsp),lapw%k3(vec,jsp)/)
299

300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322
                  tmk=tpi_const*DOT_PRODUCT(lapw%k1(:,jsp),atoms%taual(:,natom))
                  phase = CMPLX(COS(tmk),SIN(tmk))
                  fkp=MATMUL(fk,cell%bmat)
                  CALL ylm4(3,fkp,ylm)
                  DO locmm=1,mqnum
                     lmm=ll1+locmm
                     loccoffs(locmm,locm)=ci**angmom*phase*CONJG(ylm(lmm+1))
                  ENDDO
               ENDDO !locm
               !*************************************************************
               !        If locs are entangled by symmetry, the number of locs
               !        that have to be treated at the same time is larger.
               !        => l_spacetwo: two locs, which are NOT in the
               !        xy-plane and which are not related by inversion, are
               !        entangled by z-reflexion.
               !        => l_spacefour: two pairs the partners of which are
               !        connected by inversion are interrelated by z-reflection.
               !*************************************************************
               IF(l_spacetwo.OR.l_spacefour)THEN
                  DO locm=1,locnum
                     pos=basindex(zrefatom,loc)-1+locm
                     vec=kveclo(pos)
                     fk(:)=bkpt(:)+(/lapw%k1(vec,jsp),lapw%k2(vec,jsp),lapw%k3(vec,jsp)/)
323

324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444
                     tmk=tpi_const*DOT_PRODUCT(lapw%k1(:,jsp),atoms%taual(:,zrefatom))
                     phase = CMPLX(COS(tmk),SIN(tmk))
                     fkp=MATMUL(fk,cell%bmat)
                     CALL ylm4(3,fkp,ylm)
                     DO locmm=1,mqnum
                        lmm=ll1+locmm
                        loccoffs(locnum+locmm,locnum+locm)= ci**angmom*phase*CONJG(ylm(lmm+1))
                     ENDDO
                  ENDDO !locm
                  jump(zrefatom)=jump(zrefatom)+locnum
                  IF(l_spacetwo)locnum=locnum*2
               ENDIF !l_spacetwo
               !*************************************************************
               !        If locs are entangled by symmetry, loccoffs has a larger
               !        number of components for a given loc.
               !        l_planetwo => two locs in xy-plane are entangled
               !        by inversion symmetry.
               !        l_spacetwoinv => two z-reflexion-symmetric atoms
               !        are entangled also by inversion.
               !        l_crosstwo => two inversion-symmetric atoms
               !        are not directly connected by z-reflexion.
               !        l_spacefour =>  two pairs the partners of which are
               !        connected by inversion are interrelated by z-reflection.
               !*************************************************************
               IF(l_planetwo.OR.l_crosstwo.OR.l_spacetwoinv .OR.l_spacefour) THEN
                  DO locm=1,locnum
                     DO locmm=1,mqnum
                        loccoffs(mqnum+locmm,locm)=(-1)**(locmm-1)* CONJG(loccoffs(mqnum+1-locmm,locm))
                     ENDDO! locmm
                  ENDDO ! locm
               ENDIF ! l_planetwo

               !******************************************************************
               !        If four atoms are entangled, something remains to be added
               !        here:
               !******************************************************************
               IF(l_spacefour)THEN
                  DO locm=1,locnum
                     DO locmm=1,mqnum
                        loccoffs(locnum+mqnum+locmm,locnum+locm)=(-1)**(locmm-1)*&
                             CONJG(loccoffs(locnum+mqnum+1-locmm,locnum+locm))
                     ENDDO! locmm
                  ENDDO ! locm
                  locnum=locnum*2
               ENDIF
               !******************************************************************
               !        Find linear combinations of locs that are eigenfunctions
               !        of the z-reflection operation by solving linear system
               !        of equations. Write the transformation into locrec-matrix.
               !******************************************************************
               CALL CPP_LAPACK_cgetrf(locnum,locnum,loccoffs,28,ipiv,info)
               IF(info /= 0)  CALL juDFT_error("cgetrf",calledby ="locrectify")
               DO loopi=1,locnum
                  numorec=numorec+1
                  kindlocrec(numorec)=coffkind(loopi,coffindex)
                  IF(l_invsatswap)THEN
                     reccoffs(1:locnum/2)=coeffi(1:locnum/2,loopi,coffindex)
                     reccoffs(locnum/2+1:locnum/4*3)= coeffi(locnum/4*3+1:locnum,loopi,coffindex)
                     reccoffs(locnum/4*3+1:locnum)= coeffi(locnum/2+1:locnum/4*3,loopi,coffindex)
                  ELSE
                     reccoffs(1:locnum)=coeffi(1:locnum,loopi,coffindex)
                  ENDIF
                  CALL CPP_LAPACK_cgetrs('N',locnum,1,loccoffs,28,ipiv, reccoffs,28,info)
                  IF(info /= 0)  CALL juDFT_error("cgetrs", calledby="locrectify")
                  IF (l_real) THEN
                     IF(ANY(ABS(AIMAG(reccoffs(1:locnum))).GT.1e-6)) THEN
                        !In the real version of the code, the hamiltonian and overlap
                        !matrices are real. Hence the transformation matrices that make
                        !the locs become eigenfunctions of the z-reflection operation
                        !have to be real also. Check this here. 
                        PRINT*,"Sorry!"
                        PRINT*,"The transformation is not purely real."
                        PRINT*,"=> Something is wrong here."
                        CALL juDFT_error("unfortunately complex",calledby ="locrectify")
                     ENDIF
                     IF(l_spacefour)THEN
                        locrec_r(bas+1:bas+locnum/2,numorec)=REAL(reccoffs(:locnum/2))
                        DO locmm=1,locnum/2
                           locrec_r(basindex(zrefatom,loc)-1+locmm,numorec)= REAL(reccoffs(locmm+locnum/2))
                        ENDDO

                     ELSE
                        locrec_r(bas+1:bas+locnum,numorec)=REAL(reccoffs(:locnum))

                     ENDIF
                  ELSE
                     IF(l_spacetwo)THEN
                        locrec_c(bas+1:bas+locnum/2,numorec)=reccoffs(1:locnum/2)
                        locrec_c(basindex(zrefatom,loc):basindex(zrefatom,loc)+locnum/2-1,numorec)=&
                             reccoffs(1+locnum/2:locnum)
                     ELSE
                        locrec_c(bas+1:bas+locnum,numorec)=reccoffs(1:locnum)
                     ENDIF
                  ENDIF
               ENDDO !loopi
               IF(l_spacetwo) locnum=locnum/2
               IF(l_spacefour) locnum=locnum/2
               bas=bas+locnum
            ENDDO !loc
         ENDDO !natom
         nalo=naup+1
      ENDDO !ind


      !*********************************************************
      !     Determine number of even and odd locs.
      !     Prepare sort-arrays.
      !*********************************************************
      evenlocs=0
      oddlocs=0
      DO loopi=1,atoms%nlotot
         IF(kindlocrec(loopi)) THEN
            evenlocs=evenlocs+1
            evenlocindex(evenlocs)=loopi
         ELSE
            oddlocs=oddlocs+1
            oddlocindex(oddlocs)=loopi
         ENDIF
      ENDDO
    END SUBROUTINE locrectify
END MODULE m_locrectify
445