gen_wavf.F90 12.3 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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     This module generates the cmt coefficients and eigenvectors z   !
!     at all kpoints nkpt from the irreducible kpoints nkpti          !
!     and writes them out in cmt and z, respectively.                 !
!                                                 M.Betzinger(09/07)  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
14
MODULE m_gen_wavf
15

16
CONTAINS
17

18
   SUBROUTINE gen_wavf (nkpti,kpts,it,sym,atoms,el_eig,ello_eig,cell,dimension,hybrid,vr0,&
19
                        hybdat,noco,oneD,mpi,input,jsp,zmat)
20 21 22 23 24 25 26 27 28 29 30

      ! nkpti      ::     number of irreducible k-points
      ! nkpt       ::     number of all k-points 

      USE m_radfun
      USE m_radflo
      USE m_abcof
      USE m_trafo     ,ONLY: waveftrafo_genwavf
      USE m_util      ,ONLY: modulo1
      USE m_olap
      USE m_types
31
      USE m_hyb_abcrot
Daniel Wortmann's avatar
Daniel Wortmann committed
32
      USE m_io_hybrid
33

34
      IMPLICIT NONE
35

36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
      TYPE(t_hybdat),    INTENT(INOUT) :: hybdat
      TYPE(t_mpi),       INTENT(IN)    :: mpi
      TYPE(t_dimension), INTENT(IN)    :: dimension
      TYPE(t_oneD),      INTENT(IN)    :: oneD
      TYPE(t_hybrid),    INTENT(IN)    :: hybrid
      TYPE(t_input),     INTENT(IN)    :: input
      TYPE(t_noco),      INTENT(IN)    :: noco
      TYPE(t_sym),       INTENT(IN)    :: sym
      TYPE(t_cell),      INTENT(IN)    :: cell
      TYPE(t_kpts),      INTENT(IN)    :: kpts
      TYPE(t_atoms),     INTENT(IN)    :: atoms
      TYPE(t_mat),       INTENT(IN)    :: zmat(:) !for all kpoints 

      INTEGER,           INTENT(IN)    :: nkpti, it
      INTEGER,           INTENT(IN)    :: jsp

      REAL,              INTENT(IN)    :: vr0(:,:,:)!(jmtd,ntype,jspd)
      REAL,              INTENT(IN)    :: el_eig(0:atoms%lmaxd,atoms%ntype)
      REAL,              INTENT(IN)    :: ello_eig(atoms%nlod,atoms%ntype)

      ! local scalars
      INTEGER                 :: ilo,idum,m,irecl_cmt,irecl_z
      COMPLEX                 :: cdum
Daniel Wortmann's avatar
Daniel Wortmann committed
59
      TYPE(t_mat)             :: zhlp
60 61 62
      INTEGER                 :: ikpt0,ikpt,itype,iop,ispin,ieq,indx,iatom
      INTEGER                 :: i,j,l ,ll,lm,ng,ok
      COMPLEX                 :: img=(0d0,1d0)
63

64 65
      INTEGER                 :: nodem,noded
      REAL                    :: wronk
66

67 68
      INTEGER                 :: lower, upper
      LOGICAL                 :: found
69

70 71 72 73 74 75 76
      ! local arrays
      INTEGER                 :: rrot(3,3,sym%nsym)
      INTEGER                 :: map_lo(atoms%nlod)
      INTEGER                 :: iarr(0:atoms%lmaxd,atoms%ntype)
      COMPLEX,ALLOCATABLE     :: acof(:,:,:),bcof(:,:,:),ccof(:,:,:,:)
      
      COMPLEX,ALLOCATABLE     :: cmt(:,:,:),cmthlp(:,:,:)
77

78
      REAL                    :: vr(atoms%jmtd,atoms%ntype,input%jspins)
79
      REAL,ALLOCATABLE        :: f(:,:,:),df(:,:,:)
80

81 82 83
      REAL                    :: flo(atoms%jmtd,2,atoms%nlod)
      REAL                    :: uuilon(atoms%nlod,atoms%ntype),duilon(atoms%nlod,atoms%ntype)
      REAL                    :: ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)
Daniel Wortmann's avatar
Daniel Wortmann committed
84
   
85
      REAL                    :: bkpt(3)
86

87 88
!     local arrays for abcof1
!      COMPLEX                 ::  a(nvd,0:lmd,natd,nkpti),b(nvd,0:lmd,natd,nkpti)
89 90


91 92
      TYPE(t_lapw)  :: lapw(kpts%nkptf)
      TYPE(t_usdus) :: usdus
93

94
      CALL usdus%init(atoms,input%jspins)
95
      CALL zhlp%alloc(zmat(1)%l_real,zmat(1)%matsize1,zmat(1)%matsize2)
Daniel Wortmann's avatar
Daniel Wortmann committed
96 97

      
98
      ! setup rotations in reciprocal space
99 100 101 102 103 104
      DO iop = 1, sym%nsym
         IF(iop.LE.sym%nop) THEN
            rrot(:,:,iop) = transpose(sym%mrot(:,:,sym%invtab(iop)))
         ELSE
            rrot(:,:,iop) = -rrot(:,:,iop-sym%nop)
         END IF
105 106 107 108
      END DO

      ! generate G-vectors, which fulfill |k+G|<rkmax
      ! for all k-points
109 110
      DO ikpt = 1, kpts%nkptf
         CALL lapw(ikpt)%init(input,noco,kpts,atoms,sym,ikpt,cell,sym%zrfs)
111 112
      END DO

113
      ! set spherical component of the potential from the previous iteration vr
114
      vr = vr0
115 116 117 118 119 120 121 122 123 124 125 126 127


!       ALLOCATE ( z_out(nbasfcn,neigd,nkpti),stat=ok )
!       IF ( ok .ne. 0) STOP 'gen_wavf: failure allocation z'
!       z_out = 0
!       z_out(:,:,:nkpti) = z_in


      ! calculate radial basis functions belonging to the
      ! potential vr stored in bas1 and bas2
      ! bas1 denotes the large component
      ! bas2    "     "  small component

128 129 130
      ALLOCATE(f(atoms%jmtd,2,0:atoms%lmaxd),df(atoms%jmtd,2,0:atoms%lmaxd))
      f = 0
      df = 0
131
      iarr = 2
132 133 134 135
      DO itype = 1, atoms%ntype
         IF (mpi%irank == 0) WRITE (6,FMT=8000) itype
         ng = atoms%jri(itype)
         DO l=0,atoms%lmax(itype)
136 137 138
            CALL radfun(l,itype,jsp,el_eig(l,itype),vr(:,itype,jsp),atoms,f(:,:,l),df(:,:,l),usdus,nodem,noded,wronk)
            IF (mpi%irank == 0 ) WRITE (6,FMT=8010) l,el_eig(l,itype),usdus%us(l,itype,jsp),usdus%dus(l,itype,jsp),nodem,&
                                                    usdus%uds(l,itype,jsp),usdus%duds(l,itype,jsp),noded,usdus%ddn(l,itype,jsp),wronk
139 140 141 142 143 144

            hybdat%bas1(1:ng,1,l,itype) =  f(1:ng,1,l)
            hybdat%bas2(1:ng,1,l,itype) =  f(1:ng,2,l)
            hybdat%bas1(1:ng,2,l,itype) = df(1:ng,1,l)
            hybdat%bas2(1:ng,2,l,itype) = df(1:ng,2,l)

145 146 147 148
            hybdat%bas1_MT(1,l,itype)   = usdus%us(l,itype,jsp)
            hybdat%drbas1_MT(1,l,itype) = usdus%dus(l,itype,jsp)
            hybdat%bas1_MT(2,l,itype)   = usdus%uds(l,itype,jsp)
            hybdat%drbas1_MT(2,l,itype) = usdus%duds(l,itype,jsp)
149 150 151 152 153 154 155 156 157
         END DO

         IF (atoms%nlo(itype).GE.1) THEN
            CALL radflo(atoms,itype,jsp,ello_eig,vr(:,itype,jsp),f,df,mpi,usdus,uuilon,duilon,ulouilopn,flo)

            DO ilo=1,atoms%nlo(itype)
               iarr(atoms%llo(ilo,itype),itype) = iarr(atoms%llo(ilo,itype),itype) + 1
               hybdat%bas1(1:ng,iarr(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),itype) = flo(1:ng,1,ilo)
               hybdat%bas2(1:ng,iarr(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),itype) = flo(1:ng,2,ilo)
158 159
               hybdat%bas1_MT(iarr(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),itype) = usdus%ulos(ilo,itype,jsp)
               hybdat%drbas1_MT(iarr(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),itype) = usdus%dulos(ilo,itype,jsp)
160 161
            END DO
         END IF
162 163
      END DO
      DEALLOCATE (f,df)
164

165 166
#if CPP_DEBUG
      ! consistency check
167
      IF(.not.all(iarr.eq.hybrid%nindx)) STOP 'gen_wavf: counting error'
168 169
#endif

170 171 172 173
      8000 FORMAT (1x,/,/,' wavefunction parameters for atom type',i3,':',/,t32,'radial function',t79,&
                   'energy derivative',/,t3,'l',t8,'energy',t26,'value',t39,'derivative',t53,&
                   'nodes',t68,'value',t81,'derivative',t95,'nodes',t107,'norm',t119,'wronskian')
      8010 FORMAT (i3,f10.5,2 (5x,1p,2e16.7,i5),1p,2e16.7)
174 175 176 177 178

      ! determine boundaries for parallel calculations
      lower = 1
      upper = nkpti
      found = .false.
179 180
#ifdef CPP_MPI
      DO ikpt = 1, nkpti
181
         IF (.NOT.found) THEN
182 183
            lower = ikpt
            found = .true.
184 185 186 187 188 189 190
         END IF
      END DO
#else
      found = .true.
#endif
      IF (.NOT.found) THEN
         upper = 0
191 192 193 194 195 196
      END IF

      ! calculate wavefunction expansion in the the MT region
      ! (acof,bcof,ccof) and APW-basis coefficients
      ! (a,b,bascofold_lo) at irred. kpoints

197 198 199 200 201 202 203 204 205 206
      ALLOCATE(acof(dimension%neigd,0:dimension%lmd,atoms%nat),stat=ok)
      IF (ok.NE.0) STOP 'gen_wavf: failure allocation acof'
      ALLOCATE(bcof(dimension%neigd,0:dimension%lmd,atoms%nat),stat=ok)
      IF (ok.NE.0) STOP 'gen_wavf: failure allocation bcof'
      ALLOCATE(ccof(-atoms%llod:atoms%llod,dimension%neigd,atoms%nlod,atoms%nat),stat=ok)
      IF (ok.NE.0) STOP 'gen_wavf: failure allocation ccof'
      ALLOCATE (cmt(dimension%neigd,hybrid%maxlmindx,atoms%nat),stat=ok)
      IF (ok.NE.0) STOP 'gen_wavf: Failure allocation cmt'
      ALLOCATE (cmthlp(dimension%neigd,hybrid%maxlmindx,atoms%nat),stat=ok)
      IF (ok.NE.0) STOP 'gen_wavf: failure allocation cmthlp'
207 208 209

      DO ikpt0 = lower, upper

210
         acof = 0; bcof = 0; ccof = 0
211

212 213 214 215 216 217
         ! abcof calculates the wavefunction coefficients
         ! stored in acof,bcof,ccof
         lapw(ikpt0)%nmat = lapw(ikpt0)%nv(jsp) + atoms%nlotot
         CALL abcof(input,atoms,sym,cell,lapw(ikpt0),hybrid%nbands(ikpt0),usdus,noco,jsp,&!hybdat%kveclo_eig(:,ikpt0),&
                   oneD,acof(: hybrid%nbands(ikpt0),:,:),bcof(: hybrid%nbands(ikpt0),:,:),&
                   ccof(:,: hybrid%nbands(ikpt0),:,:),zmat(ikpt0))
218 219 220 221 222 223 224 225 226 227 228

! call was ...
          ! gpt(1,:,:,ikpt0),gpt(2,:,:,ikpt0),&
          ! gpt(3,:,:,ikpt0),ngpt(:,ikpt0),&!k1hlp,k2hlp,k3hlp,nvhlp,&
          !    ngpt(jsp,ikpt0)+nbands(ikpt0),z(:,:,ikpt0),&!nvhlp(jsp)+ &
          !   &usdus,&
          !    noco,&
          !    jsp,kveclo_eig(:ikpt0),oneD,oneD,&
          !    acof(:nbands(ikpt0),:,:),&
          !    bcof(:nbands(ikpt0),:,:),ccof(:,:nbands(ikpt0),:,:) )

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
         ! MT wavefunction coefficients are calculated in a local coordinate system rotate them in the global one

         CALL hyb_abcrot(hybrid,atoms,hybrid%nbands(ikpt0),sym,cell,oneD,acof(: hybrid%nbands(ikpt0),:,:),&
                         bcof(: hybrid%nbands(ikpt0),:,:),ccof(:,: hybrid%nbands(ikpt0),:,:))

         ! decorate acof, bcof, ccof with coefficient i**l and store them
         ! in the field cmt(neigd,nkpt,maxlmindx,nat), i.e.
         ! where maxlmindx subsumes l,m and nindx

         cmt = 0
         iatom = 0
         DO itype = 1, atoms%ntype
            DO ieq = 1, atoms%neq(itype)
               iatom = iatom + 1
               indx = 0
               DO l = 0, atoms%lmax(itype)
                  ll = l*(l+1)
                  cdum = img**l

                  ! determine number of local orbitals with quantum number l
                  ! map returns the number of the local orbital of quantum
                  ! number l in the list of all local orbitals of the atom type
                  idum   = 0
                  map_lo = 0
                  IF (hybrid%nindx(l,itype).GT.2) THEN
                     DO j = 1, atoms%nlo(itype)
                        IF (atoms%llo(j,itype).EQ.l) THEN
                           idum = idum + 1
                           map_lo(idum) = j
                        END IF
                     END DO
260 261
                  END IF

262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278
                  DO M = -l, l
                     lm = ll + M
                     DO i = 1, hybrid%nindx(l,itype)
                        indx = indx + 1
                        IF(i.EQ.1) THEN
                           cmt(:,indx,iatom) = cdum * acof(:,lm,iatom)
                        ELSE IF (i.EQ.2) THEN
                           cmt(:,indx,iatom) = cdum * bcof(:,lm,iatom)
                        ELSE
                           idum = i - 2
                           cmt(:,indx,iatom) = cdum * ccof(M,:,map_lo(idum),iatom)
                        END IF
                     END DO
                  END DO
               END DO
            END DO
         END DO
279

280 281
         ! write cmt at irreducible k-points in direct-access file cmt
         CALL write_cmt(cmt,ikpt0)
282
         CALL zhlp%alloc(zmat(ikpt0)%l_real,zmat(ikpt0)%matsize1,zmat(ikpt0)%matsize2)
Daniel Wortmann's avatar
Daniel Wortmann committed
283
        
284 285 286 287 288 289 290 291 292
         IF (zhlp%l_real) THEN
            zhlp%data_r = zmat(ikpt0)%data_r
         ELSE
            zhlp%data_c = zmat(ikpt0)%data_c
         END IF
         CALL write_z(zhlp,ikpt0)

         ! generate wavefunctions coefficients at all k-points from
         ! irreducible k-points
293

294 295 296 297 298 299 300 301 302 303 304
         DO ikpt = 1, kpts%nkptf
            IF ((kpts%bkp(ikpt).EQ.ikpt0).AND.(ikpt0.NE.ikpt)) THEN
               iop = kpts%bksym(ikpt)
               CALL waveftrafo_genwavf(cmthlp,zhlp%data_r,zhlp%data_c,cmt(:,:,:),zmat(1)%l_real,zmat(ikpt0)%data_r(:,:),&
                                       zmat(ikpt0)%data_c(:,:),ikpt0,iop,atoms,hybrid,kpts,sym,jsp,dimension,&
                                       hybrid%nbands(ikpt0),cell,lapw(ikpt0),lapw(ikpt),.true.)

               CALL write_cmt(cmthlp,ikpt)
               CALL write_z(zhlp,ikpt)
            END IF
         END DO  !ikpt
305 306
      END DO !ikpt0

307 308
      DEALLOCATE(acof,bcof,ccof)
      DEALLOCATE(cmt,cmthlp)
309

310
   END SUBROUTINE gen_wavf
311

312
END MODULE m_gen_wavf