gen_wavf.F90 15.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     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)  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
      MODULE m_gen_wavf

      CONTAINS

      SUBROUTINE gen_wavf (&
     &          nkpti,kpts,it,sym,&
     &          atoms,el_eig,ello_eig,cell,&
     &          dimension,hybrid,vr0,&
Daniel Wortmann's avatar
Daniel Wortmann committed
16
     &          hybdat,&
17
     &          noco,oneD,mpi,irank2,&
Daniel Wortmann's avatar
Daniel Wortmann committed
18
19
     &          input,jsp,&
     &          zmat)
20
21
22
23
24
25
26
27
28
29
30
31
32
33


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

      USE m_apws
      USE m_radfun
      USE m_radflo
      USE m_abcof
      USE m_trafo     ,ONLY: waveftrafo_genwavf
      USE m_util      ,ONLY: modulo1
      USE m_setabc1locdn
      USE m_olap
      USE m_types
Daniel Wortmann's avatar
Daniel Wortmann committed
34
      USE m_abcrot
35
      IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
36
37
38

      TYPE(t_hybdat),INTENT(INOUT)   :: hybdat
      TYPE(t_mpi),INTENT(IN)         :: mpi
39
      TYPE(t_dimension),INTENT(IN)   :: dimension
Daniel Wortmann's avatar
Daniel Wortmann committed
40
41
42
43
44
45
46
47
      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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

!     - - scalars - -
      INTEGER,INTENT(IN)      ::  nkpti ,it    

!     scalars for abcoff

!     scalars for apws
      INTEGER,INTENT(IN)      :: jsp

!     - - arrays - -
      INTEGER,INTENT(IN)      ::  irank2(nkpti)


      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)
      TYPE(t_zmat),INTENT(IN) :: zmat(:) !for all kpoints 

  !     - - local scalars - - 
Daniel Wortmann's avatar
Daniel Wortmann committed
68
      INTEGER                 ::  ilo,idum ,m,maxlmindx
69
70
71
72
73
74
75
76
77
78
79
80
81
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
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
      REAL                    ::  rdum,merror
      COMPLEX                 ::  cdum,cdum1,cdum2

!     local scalars for apws
      INTEGER                 ::  nred
      INTEGER                 ::  ikpt0,ikpt,ikpt1,iband,itype,iop,&
     &                            ispin,ieq,ic,indx,iatom
      INTEGER                 ::  i,j,l ,ll,lm,nrkpt,ng,ok
      REAL                    ::  time1,time2,time3,time4
      COMPLEX                 ::  img=(0d0,1d0)

!     local scalars for radfun
      INTEGER                 ::  nodem,noded
      REAL                    ::  wronk

!     reduced dimension for parallel calculation
      INTEGER                 ::  lower, upper
      LOGICAL                 ::  found

!     - - local arrays - -
      INTEGER                 ::  rrot(3,3,sym%nsym)
      INTEGER                 ::  map_lo(atoms%nlod)
      INTEGER                 ::  iarr(0:atoms%lmaxd,atoms%ntype)
      REAL                    ::  rotkpt(3)
      COMPLEX,ALLOCATABLE     ::  acof(:,:,:),bcof(:,:,:),ccof(:,:,:,:)
      REAL,ALLOCATABLE        ::  zhlp_r(:,:)
      COMPLEX,ALLOCATABLE     ::  zhlp_c(:,:)

      COMPLEX,ALLOCATABLE     ::  cmt(:,:,:),cmthlp(:,:,:)


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

!     local arrays for radfun
      REAL                    ::  vr(atoms%jmtd,atoms%ntype,dimension%jspd)
      REAL,ALLOCATABLE        ::  f(:,:,:),df(:,:,:)


!     local arrays for radflo
      REAL                    ::  flo(atoms%jmtd,2,atoms%nlod)
      REAL                    ::  uloulopn(atoms%nlod,atoms%nlod,atoms%ntype)
      REAL                    ::  uuilon(atoms%nlod,atoms%ntype),duilon(atoms%nlod,atoms%ntype)
      REAL                    ::  ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)

!     local arrays for setabc1locdn
      INTEGER                 ::  nbasf0(atoms%nlod,atoms%nat),nkvec(atoms%nlod,atoms%nat)
      INTEGER                 ::  kvec(2*(2*atoms%llod+1)  )
      LOGICAL                 ::  enough(atoms%nat)


!     local arrays for apws
      INTEGER                 ::  matind(dimension%nbasfcn,2)
      INTEGER                 ::  gpthlp1(3,dimension%nvd,dimension%jspd),nvhlp1(dimension%jspd)
      INTEGER                 ::  gpthlp2(3,dimension%nvd,dimension%jspd),nvhlp2(dimension%jspd)
      INTEGER                 ::  k1hlp(dimension%nvd,dimension%jspd),k2hlp(dimension%nvd,dimension%jspd),&
     &                            k3hlp(dimension%nvd,dimension%jspd),nvhlp(dimension%jspd)
      REAL                    :: bkpt(3)

      REAL,ALLOCATABLE        ::  olapmt(:,:,:,:)
#if ( defined(CPP_INVERSION) )
      REAL,ALLOCATABLE        ::  olappw(:,:)
#else
      COMPLEX,ALLOCATABLE     ::  olappw(:,:)
#endif
      INTEGER                 ::  irecl_cmt,irecl_z

      INTEGER                 ::  gpt(3,dimension%nvd,dimension%jspd,kpts%nkpt),ngpt(dimension%jspd,kpts%nkpt)

      TYPE(t_lapw) :: lapw
      TYPE(t_usdus):: usdus

Daniel Wortmann's avatar
Daniel Wortmann committed
141
142
143
      CALL CPU_TIME(time1)
      call usdus%init(atoms,dimension%jspd)
     
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
      ! setup rotations in reciprocal space
      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
      END DO

      ! generate G-vectors, which fulfill |k+G|<rkmax
      ! for all k-points
      DO ikpt=1,kpts%nkpt
        CALL apws(dimension,input,noco,&
     &            kpts,ikpt,cell,sym%zrfs,&
     &            1,jsp,bkpt,lapw,matind,nred)

      END DO

      ! read in spherical component of the potential from
      ! the previous iteration
      ! eigen_hf writes spherical component in file vr0
      ! if it = 1 vr0 is identical to the potential of the
      ! previous calculation

      IF( it .ne. 1) THEN
        OPEN(unit=220,file='vr0',form='unformatted')
        DO ispin=1,dimension%jspd
          DO itype=1,atoms%ntype
            DO i=1,atoms%jmtd
              READ(220) vr(i,itype,ispin)
            END DO
          END DO
        END DO
        CLOSE(220)
      ELSE
        vr = vr0
      END IF


!       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

      ALLOCATE( f(atoms%jmtd,2,0:atoms%lmaxd), df(atoms%jmtd,2,0:atoms%lmaxd) )
      f    = 0
      df   = 0
      iarr = 2
      DO itype=1,atoms%ntype
        if ( mpi%irank == 0 ) WRITE (6,FMT=8000) itype
        ng = atoms%jri(itype)
        DO l=0,atoms%lmax(itype)
           CALL radfun(l,itype,1,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,1),usdus%dus(l,itype,1),nodem,usdus%uds(l,itype,1),&
               usdus%duds(l,itype,1),noded,usdus%ddn(l,itype,1),wronk

Daniel Wortmann's avatar
Daniel Wortmann committed
208
209
210
211
          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)
212

Daniel Wortmann's avatar
Daniel Wortmann committed
213
214
215
216
            hybdat%bas1_MT(1,l,itype) =   usdus%us(l,itype,1)
          hybdat%drbas1_MT(1,l,itype) =  usdus%dus(l,itype,1)
            hybdat%bas1_MT(2,l,itype) =  usdus%uds(l,itype,1)
          hybdat%drbas1_MT(2,l,itype) = usdus%duds(l,itype,1)
217
218
219
220
221
222
223
224
225
226
227
        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
Daniel Wortmann's avatar
Daniel Wortmann committed
228
            hybdat%bas1(1:ng,iarr(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),itype)&
229
     &                                                = flo(1:ng,1,ilo)
Daniel Wortmann's avatar
Daniel Wortmann committed
230
            hybdat%bas2(1:ng,iarr(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),itype)&
231
232
     &                                                = flo(1:ng,2,ilo)

Daniel Wortmann's avatar
Daniel Wortmann committed
233
              hybdat%bas1_MT(iarr(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),itype)&
234
     &                                                = usdus%ulos(ilo,itype,1)
Daniel Wortmann's avatar
Daniel Wortmann committed
235
            hybdat%drbas1_MT(iarr(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),itype)&
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
283
284
285
286
287
288
289
290
291
292
293
     &                                                = usdus%dulos(ilo,itype,1)
          END DO

        END IF
      END DO
      DEALLOCATE (f,df)
#if CPP_DEBUG
      ! consistency check
      IF( .not. all(iarr .eq. hybrid%nindx ) ) STOP 'gen_wavf: counting error'
#endif

 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)


      ! determine boundaries for parallel calculations
      lower = 1
      upper = nkpti
      found = .false.
#     ifdef CPP_MPI
        DO ikpt = 1, nkpti
          IF ( irank2(ikpt) == 0 .AND. .not. found ) THEN
            lower = ikpt
            found = .true.
          ELSE IF ( irank2(ikpt) /= 0 .AND. found ) THEN
            upper = ikpt-1
            EXIT
          END IF
        END DO
#     else
        found = .true.
#     endif
      IF ( .not. found ) THEN
        upper = 0
      END IF


      maxlmindx = maxval(&
     &       (/ ( sum( (/ (hybrid%nindx(l,itype)*(2*l+1),l=0,atoms%lmax(itype)) /) ),&
     &                                              itype=1,atoms%ntype) /) )

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

!       CALL cpu_time(time2)
!       WRITE(*,*) 'time for generating radial functions',time2-time1

      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'
Daniel Wortmann's avatar
Daniel Wortmann committed
294
      ALLOCATE ( cmt(dimension%neigd,maxlmindx,atoms%nat), stat=ok)
295
      IF(  ok .ne. 0 ) STOP 'gen_wavf: Failure allocation cmt'
Daniel Wortmann's avatar
Daniel Wortmann committed
296
      ALLOCATE ( cmthlp(dimension%neigd,maxlmindx,atoms%nat), stat=ok)
297
298
299
300
301
302
303
304
      IF( ok .ne. 0) STOP 'gen_wavf: failure allocation cmthlp'
      if (zmat(1)%l_real) THEN
         ALLOCATE ( zhlp_r(dimension%nbasfcn,dimension%neigd), stat=ok)
      ELSE
         ALLOCATE ( zhlp_c(dimension%nbasfcn,dimension%neigd), stat=ok)
      ENDIF
      IF( ok .ne. 0) STOP 'gen_wavf: failure allocation zhlp'

Daniel Wortmann's avatar
Daniel Wortmann committed
305
      irecl_cmt = dimension%neigd*maxlmindx*atoms%nat*16
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
      OPEN(unit=777,file='cmt',form='unformatted',access='direct',&
     &     recl=irecl_cmt)

#     ifdef CPP_INVERSION
        irecl_z   =  dimension%nbasfcn*dimension%neigd*8
#     else
        irecl_z   =  dimension%nbasfcn*dimension%neigd*16
#     endif
      OPEN(unit=778,file='z',form='unformatted',access='direct',&
     &     recl=irecl_z)

      DO ikpt0 = lower, upper

        acof = 0; bcof = 0; ccof = 0

        ! abcof calculates the wavefunction coefficients
        ! stored in acof,bcof,ccof
        CALL abcof(&
Daniel Wortmann's avatar
Daniel Wortmann committed
324
325
326
              input,atoms,hybdat%nbands(ikpt0),sym, cell, Kpts%bk(:,ikpt0), lapw, &
              ngpt(jsp,ikpt0)+hybdat%nbands(ikpt0),usdus,noco,jsp,hybdat%kveclo_eig(:,ikpt0),&
              oneD,acof(: hybdat%nbands(ikpt0),:,:),bcof(: hybdat%nbands(ikpt0),:,:),ccof(:,: hybdat%nbands(ikpt0),:,:),&
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
              zmat(ikpt))
        

! 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),:,:) )

        ! MT wavefunction coefficients are calculated in a local coordinate system
        ! rotate them in the global one

        CALL abcrot(&
Daniel Wortmann's avatar
Daniel Wortmann committed
344
345
                atoms,hybdat%nbands(ikpt0),&
                 sym,&
346
                cell,oneD,&
Daniel Wortmann's avatar
Daniel Wortmann committed
347
348
                acof(: hybdat%nbands(ikpt0),:,:),bcof(: hybdat%nbands(ikpt0),:,:),&
                ccof(:,: hybdat%nbands(ikpt0),:,:) )
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

!       CALL cpu_time(time3)

!       WRITE(*,*) 'time for abcoff and abcrot',time3-time2


!       CALL cpu_time(time2)


        ! 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
              END IF



              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)
                  ELSEIF ( 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

!       CALL cpu_time(time3)

        ! write cmt at irreducible k-points in direct-access file cmt
        WRITE(777,rec=ikpt0) cmt

        ! write z at irreducible k-points in direct-access file z
        IF (zmat(1)%l_real) THEN
           WRITE(778,rec=ikpt0) zmat(ikpt0)%z_r(:,:)
        ELSE
           WRITE(778,rec=ikpt0) zmat(ikpt0)%z_c(:,:)
        ENDIF
       

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

        DO ikpt=1,kpts%nkpt
          IF ( kpts%bkp(ikpt) .eq. ikpt0 .and. ikpt0 .ne. ikpt ) THEN
            iop = kpts%bksym(ikpt)


 
            
              CALL waveftrafo_genwavf( cmthlp,zhlp_r,zhlp_c,&
     &                 cmt(:,:,:),zmat(1)%l_real,zmat(ikpt0)%z_r(:,:),zmat(ikpt0)%z_c(:,:),ikpt0,iop,atoms,&
     &                 hybrid,kpts,sym,&
Daniel Wortmann's avatar
Daniel Wortmann committed
433
     &                 jsp,dimension,hybdat%nbands(ikpt0),&
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
     &                 cell,gpt(:,:ngpt(jsp,ikpt0),jsp,ikpt0),&
     &                 ngpt(:,ikpt0),gpt(:,:ngpt(jsp,ikpt),jsp,ikpt),&
     &                 ngpt(:,ikpt),.true.)

              WRITE(777,rec=ikpt) cmthlp
              IF (zmat(1)%l_real) THEN
                 WRITE(778,rec=ikpt) zhlp_r
              ELSE
                 WRITE(778,rec=ikpt) zhlp_c
              ENDIF
          END IF
        END DO  !ikpt
      END DO !ikpt0

      DEALLOCATE( acof,bcof,ccof )
      DEALLOCATE( cmt,cmthlp)

      !close file cmt and z
      CLOSE(777)
      CLOSE(778)

      END SUBROUTINE gen_wavf

      END MODULE m_gen_wavf