hf_setup.F90 12.6 KB
Newer Older
1

Daniel Wortmann's avatar
Daniel Wortmann committed
2
MODULE m_hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
3
CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
4
  SUBROUTINE hf_setup(hybrid,input,sym,kpts,DIMENSION,atoms,mpi,noco,cell,oneD,results,jsp,eig_id_hf,&
Daniel Wortmann's avatar
Daniel Wortmann committed
5
       hybdat,irank2,it,l_real,vr0,eig_irr) 
Daniel Wortmann's avatar
Daniel Wortmann committed
6 7 8 9 10 11 12 13 14
    USE m_types
    USE m_eig66_io
    USE m_util
    USE m_apws
    USE m_checkolap
    USE m_read_core
    USE m_gen_wavf
    IMPLICIT NONE
    TYPE(t_hybrid),INTENT(INOUT):: hybrid
15
    TYPE(t_kpts),INTENT(IN)     :: kpts
Daniel Wortmann's avatar
Daniel Wortmann committed
16
    TYPE(t_dimension),INTENT(IN):: DIMENSION
17 18 19 20 21 22 23
    TYPE(t_atoms),INTENT(IN)    :: atoms
    TYPE(t_mpi),INTENT(IN)      :: mpi
    TYPE(t_noco),INTENT(IN)     :: noco
    TYPE(t_cell),INTENT(IN)     :: cell
    TYPE(t_oneD),INTENT(IN)     :: oneD
    TYPE(t_input),INTENT(IN)    :: input
    TYPE(t_sym),INTENT(IN)      :: sym
Daniel Wortmann's avatar
Daniel Wortmann committed
24
    TYPE(t_results),INTENT(INOUT):: results
25
    INTEGER,INTENT(IN)          :: irank2(:),it
Daniel Wortmann's avatar
Daniel Wortmann committed
26

27 28
    INTEGER,INTENT(IN)          :: jsp,eig_id_hf
    REAL, INTENT(IN)            :: vr0(:,:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
29 30
    LOGICAL,INTENT(IN)          :: l_real
    TYPE(t_hybdat),INTENT(INOUT):: hybdat
Daniel Wortmann's avatar
Daniel Wortmann committed
31 32
    REAL,ALLOCATABLE,INTENT(OUT):: eig_irr(:,:)
 
33

Daniel Wortmann's avatar
Daniel Wortmann committed
34
    INTEGER:: ok,nk,nrec1,i,j,ll,l1,l2,ng,itype,n,l,n1,n2,nn
Daniel Wortmann's avatar
Daniel Wortmann committed
35 36 37 38 39 40 41 42


    TYPE(t_zmat),ALLOCATABLE :: zmat(:)
    REAL,    ALLOCATABLE    ::  basprod(:)
    REAL                    ::  el_eig(0:atoms%lmaxd,atoms%ntype), ello_eig(atoms%nlod,atoms%ntype),bk(3)
    INTEGER                 ::  degenerat(DIMENSION%neigd2+1,kpts%nkpt)
    INTEGER                 :: matind(DIMENSION%nbasfcn,2),nred
    TYPE(t_lapw)            :: lapw
Daniel Wortmann's avatar
Daniel Wortmann committed
43
  
Daniel Wortmann's avatar
Daniel Wortmann committed
44 45 46 47
    LOGICAL:: skip_kpt(kpts%nkpt)
    REAL   :: g(3)
    skip_kpt=.FALSE.

48
    IF( hybrid%l_calhf ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
49
       !
50
       !             preparations for HF and hybrid functional calculation
Daniel Wortmann's avatar
Daniel Wortmann committed
51
       !
52
       CALL timestart("gen_bz and gen_wavf")
Daniel Wortmann's avatar
Daniel Wortmann committed
53

Daniel Wortmann's avatar
Daniel Wortmann committed
54
       ALLOCATE(zmat(kpts%nkptf),stat=ok)
Daniel Wortmann's avatar
Daniel Wortmann committed
55 56
       IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation z_c'
       ALLOCATE ( eig_irr(DIMENSION%neigd2,kpts%nkpt)      ,stat=ok )
Daniel Wortmann's avatar
Daniel Wortmann committed
57
       IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation eig_irr'
58
       ALLOCATE ( hybdat%kveclo_eig(atoms%nlotot,kpts%nkpt)  ,stat=ok )
Daniel Wortmann's avatar
Daniel Wortmann committed
59
       IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation hybdat%kveclo_eig'
60 61
       eig_irr = 0 ; hybdat%kveclo_eig = 0

Daniel Wortmann's avatar
Daniel Wortmann committed
62 63 64

       zmat(:)%l_real=l_real
      ! Reading the eig file
Daniel Wortmann's avatar
Daniel Wortmann committed
65
       DO nk = 1,kpts%nkpt
66
#               ifdef CPP_MPI
Daniel Wortmann's avatar
Daniel Wortmann committed
67 68
          ! jump to next k-point if this process is not present in communicator
          IF ( skip_kpt(nk) ) CYCLE
69
#               endif
Daniel Wortmann's avatar
Daniel Wortmann committed
70 71 72 73 74
          nrec1 = kpts%nkpt*(jsp-1) + nk
          zmat(nk)%nbasfcn=dimension%nbasfcn
          zmat(nk)%nbands=dimension%neigd2
          if (l_real) THEN
             ALLOCATE(zmat(nk)%z_r(dimension%nbasfcn,dimension%neigd2))
Daniel Wortmann's avatar
Daniel Wortmann committed
75
             ALLOCATE(zmat(nk)%z_c(0,0))
Daniel Wortmann's avatar
Daniel Wortmann committed
76 77
          else
             ALLOCATE(zmat(nk)%z_c(dimension%nbasfcn,dimension%neigd2))
Daniel Wortmann's avatar
Daniel Wortmann committed
78
             ALLOCATE(zmat(nk)%z_r(0,0))
Daniel Wortmann's avatar
Daniel Wortmann committed
79
          endif
Daniel Wortmann's avatar
Daniel Wortmann committed
80
          CALL read_eig(eig_id_hf,nk,jsp,el=el_eig,ello=ello_eig, neig=hybrid%ne_eig(nk),eig=eig_irr(:,nk), w_iks=results%w_iks(:,nk,jsp),kveclo=hybdat%kveclo_eig(:,nk),zmat=zmat(nk))
81
     
Daniel Wortmann's avatar
Daniel Wortmann committed
82
       END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
83 84 85 86 87 88 89 90 91 92
       !Allocate further space
       DO nk=kpts%nkpt+1,kpts%nkptf
          zmat(nk)%nbasfcn=dimension%nbasfcn
          zmat(nk)%nbands=dimension%neigd2
          if (l_real) THEN
             ALLOCATE(zmat(nk)%z_r(dimension%nbasfcn,dimension%neigd2))
          else
             ALLOCATE(zmat(nk)%z_c(dimension%nbasfcn,dimension%neigd2))
          endif
       Enddo
Daniel Wortmann's avatar
Daniel Wortmann committed
93 94 95 96 97 98 99 100 101 102 103 104 105
       !
       !determine degenerate states at each k-point
       !
       ! degenerat(i) =1  band i  is not degenerat ,
       ! degenerat(i) =j  band i  has j-1 degenart states ( i, i+1, ..., i+j)
       ! degenerat(i) =0  band i  is  degenerat, but is not the lowest band
       !                  of the group of degenerate states
       IF ( mpi%irank == 0 ) THEN
          WRITE(6,*)
          WRITE(6,'(A)') "   k-point      |   number of occupied&
               &bands  |   maximal number of bands"
       END IF
       degenerat = 1
Daniel Wortmann's avatar
Daniel Wortmann committed
106
       hybrid%nobd      = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
107
       DO nk=1 ,kpts%nkpt
108
#               ifdef CPP_MPI
Daniel Wortmann's avatar
Daniel Wortmann committed
109 110
          ! jump to next k-point if this k-point is not treated at this process
          IF ( skip_kpt(nk) ) CYCLE
111
#               endif
Daniel Wortmann's avatar
Daniel Wortmann committed
112 113
          DO i=1,hybrid%ne_eig(nk)
             DO j=i+1,hybrid%ne_eig(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
114 115
                IF( ABS(eig_irr(i,nk)-eig_irr(j,nk)) < 1E-07) THEN !0.015
                   degenerat(i,nk) = degenerat(i,nk) + 1
116
                END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
117 118 119
             END DO
          END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
120
          DO i=1,hybrid%ne_eig(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
121 122 123 124 125 126 127
             IF( degenerat(i,nk) .NE. 1 .OR. degenerat(i,nk) .NE. 0 ) &
                  degenerat(i+1:i+degenerat(i,nk)-1,nk) = 0
          END DO


          ! set the size of the exchange matrix in the space of the wavefunctions

Daniel Wortmann's avatar
Daniel Wortmann committed
128 129
          hybrid%nbands(nk)=hybrid%bands1
          IF(hybrid%nbands(nk).GT.hybrid%ne_eig(nk)) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
130
             IF ( mpi%irank == 0 ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
131
                WRITE(*,*) ' maximum for hybrid%nbands is', hybrid%ne_eig(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
132
                WRITE(*,*) ' increase energy window to obtain enough eigenvalues'
Daniel Wortmann's avatar
Daniel Wortmann committed
133
                WRITE(*,*) ' set hybrid%nbands equal to hybrid%ne_eig'
Daniel Wortmann's avatar
Daniel Wortmann committed
134
             END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
135
             hybrid%nbands(nk)=hybrid%ne_eig(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
136
          END IF
137

Daniel Wortmann's avatar
Daniel Wortmann committed
138 139 140
          DO i = hybrid%nbands(nk)-1,1,-1
             IF( (degenerat(i,nk) .GE. 1) .AND. (degenerat(i,nk)+i-1 .NE. hybrid%nbands(nk) ) ) THEN
                hybrid%nbands(nk) = i + degenerat(i,nk) - 1
Daniel Wortmann's avatar
Daniel Wortmann committed
141 142 143
                EXIT
             END IF
          END DO
144

Daniel Wortmann's avatar
Daniel Wortmann committed
145 146
          DO i = 1,hybrid%ne_eig(nk)
             IF(results%w_iks(i,nk,jsp) .GT. 0d0 ) hybrid%nobd(nk) = hybrid%nobd(nk) + 1
147

Daniel Wortmann's avatar
Daniel Wortmann committed
148
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
149
          IF (hybrid%nobd(nk)>hybrid%nbands(nk)) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
150
             CALL judft_warn("More occupied bands than total no of bands!?")
Daniel Wortmann's avatar
Daniel Wortmann committed
151
             hybrid%nbands(nk)=hybrid%nobd(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
152
          ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
153
          PRINT *,"bands:",nk, hybrid%nobd(nk),hybrid%nbands(nk),hybrid%ne_eig(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
154
       END DO
155 156

#             ifdef CPP_MPI
Daniel Wortmann's avatar
Daniel Wortmann committed
157 158 159 160 161
       ! send results for occupied bands to all processes
       sndreqd = 0 ; rcvreqd = 0
       DO nk = 1,kpts%nkpt
          IF ( skip_kpt(nk) ) THEN
             rcvreqd = rcvreqd + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
162
             CALL MPI_IRECV(hybrid%nobd(nk),1,MPI_INTEGER4, MPI_ANY_SOURCE,TAG_SNDRCV_HYBDAT%NOBD+nk, mpi,rcvreq(rcvreqd),ierr(1))
Daniel Wortmann's avatar
Daniel Wortmann committed
163 164 165 166
          ELSE
             i = MOD( mpi%irank + isize2(nk), mpi%isize )
             DO WHILE ( i < mpi%irank-irank2(nk) .OR. i >= mpi%irank-irank2(nk)+isize2(nk) )
                sndreqd = sndreqd + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
167
                CALL MPI_ISSEND(hybrid%nobd(nk),1,MPI_INTEGER4,i, TAG_SNDRCV_HYBDAT%NOBD+nk,mpi, sndreq(sndreqd),ierr(1) )
Daniel Wortmann's avatar
Daniel Wortmann committed
168 169 170 171 172 173
                i = MOD( i + isize2(nk), mpi%isize )
             END DO
          END IF
       END DO
       CALL MPI_WAITALL( rcvreqd, rcvreq, MPI_STATUSES_IGNORE, ierr(1) )
       ! Necessary to avoid compiler optimization
Daniel Wortmann's avatar
Daniel Wortmann committed
174 175
       ! Compiler does not know that hybrid%nobd is modified in mpi_waitall
       CALL MPI_GET_ADDRESS( hybrid%nobd, addr, ierr(1) )
Daniel Wortmann's avatar
Daniel Wortmann committed
176
       rcvreqd = 0
177 178 179

#             endif

Daniel Wortmann's avatar
Daniel Wortmann committed
180
       ! spread hybrid%nobd from IBZ to whole BZ
Daniel Wortmann's avatar
Daniel Wortmann committed
181 182
       DO nk = 1,kpts%nkptf
          i       = kpts%bkp(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
183
          hybrid%nobd(nk)= hybrid%nobd(i)
Daniel Wortmann's avatar
Daniel Wortmann committed
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203
       END DO

       !
       ! generate eigenvectors z and MT coefficients from the previous iteration
       ! at all k-points
       !
       CALL gen_wavf(&
            &                 kpts%nkpt,kpts,it,sym,&
            &                 atoms,el_eig,&
            &                 ello_eig,cell,DIMENSION,&
            &                 hybrid,vr0,&
            &                 hybdat,&
            &                 noco,oneD,mpi,irank2,&
            &                 input,jsp,&
            &                 zmat)

       ! generate core wave functions (-> core1/2(jmtd,hybdat%nindxc,0:lmaxc,ntype) )
       CALL corewf(atoms,jsp,input,DIMENSION,vr0,&
            &                    hybdat%lmaxcd,hybdat%maxindxc,mpi,&
            &                    hybdat%lmaxc,hybdat%nindxc,hybdat%core1,hybdat%core2,hybdat%eig_c)
204 205

#             ifdef CPP_MPI
Daniel Wortmann's avatar
Daniel Wortmann committed
206 207
       ! wait until all files are written in gen_wavf
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
208 209
#             endif

Daniel Wortmann's avatar
Daniel Wortmann committed
210 211 212 213 214 215 216 217 218 219 220 221 222 223
       !
       ! check olap between core-basis/core-valence/basis-basis
       !
       CALL checkolap(atoms,hybdat,hybrid,&
            &                      kpts%nkpt,kpts,&
            &                      DIMENSION,mpi,irank2,skip_kpt,&
            &                      input,sym,&
            &                      noco,cell,lapw,jsp)

       !
       ! set up pointer pntgpt
       !

       ! setup dimension of pntgpt
Daniel Wortmann's avatar
Daniel Wortmann committed
224
       ALLOCATE(hybdat%pntgptd(3))
Daniel Wortmann's avatar
Daniel Wortmann committed
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
       hybdat%pntgptd = 0
       DO nk = 1,kpts%nkptf
          CALL apws(DIMENSION,input,noco, kpts,nk,cell,sym%zrfs,&
               &                    1,jsp, bk,lapw,matind,nred)
          hybdat%pntgptd(1) = MAXVAL( (/ ( ABS(lapw%k1(i,jsp)),i=1,lapw%nv(jsp)), hybdat%pntgptd(1) /) )
          hybdat%pntgptd(2) = MAXVAL( (/ ( ABS(lapw%k2(i,jsp)),i=1,lapw%nv(jsp)), hybdat%pntgptd(2) /) )
          hybdat%pntgptd(3) = MAXVAL( (/ ( ABS(lapw%k3(i,jsp)),i=1,lapw%nv(jsp)), hybdat%pntgptd(3) /) )
       END DO

       ALLOCATE( hybdat%pntgpt(-hybdat%pntgptd(1):hybdat%pntgptd(1), -hybdat%pntgptd(2):hybdat%pntgptd(2),&
            &                         -hybdat%pntgptd(3):hybdat%pntgptd(3),kpts%nkptf),stat=ok )
       IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation pntgpt'
       hybdat%pntgpt = 0
       DO nk = 1,kpts%nkptf
          CALL apws( DIMENSION,input,noco, kpts,nk,cell,sym%zrfs,&
               &                     1,jsp, bk,lapw,matind,nred)
          DO i = 1,lapw%nv(jsp)
             g = (/ lapw%k1(i,jsp),lapw%k2(i,jsp),lapw%k3(i,jsp) /)
             hybdat%pntgpt(g(1),g(2),g(3),nk) = i
          END DO
       END DO

       ALLOCATE ( basprod(atoms%jmtd),stat=ok )
       IF( ok .NE. 0 )STOP 'eigen_hf: failure allocation basprod'
       ALLOCATE ( hybdat%prodm(hybrid%maxindxm1,hybrid%maxindxp1,0:hybrid%maxlcutm1,atoms%ntype), stat= ok )
       IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation hybdat%prodm'
       ALLOCATE ( hybdat%prod(hybrid%maxindxp1,0:hybrid%maxlcutm1,atoms%ntype),stat= ok )
       IF( ok .NE. 0 ) STOP 'eigen_hf: failure allocation hybdat%prod'
       basprod = 0 ; hybdat%prodm = 0 ; hybdat%prod%l1 = 0 ; hybdat%prod%l2 = 0
       hybdat%prod%n1 = 0 ; hybdat%prod%n2 = 0
255 256
       ALLOCATE(hybdat%nindxp1(0:hybrid%maxlcutm1,atoms%ntype))
       hybdat%nindxp1 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
257 258 259 260 261 262 263 264 265 266 267 268 269 270 271
       DO itype = 1,atoms%ntype
          ng = atoms%jri(itype)
          DO l2 = 0,MIN(atoms%lmax(itype),hybrid%lcutwf(itype))
             ll = l2
             DO l1 = 0,ll
                IF(ABS(l1-l2).LE.hybrid%lcutm1(itype)) THEN
                   DO n2 = 1,hybrid%nindx(l2,itype)
                      nn = hybrid%nindx(l1,itype)
                      IF(l1.EQ.l2) nn = n2
                      DO n1 = 1,nn
                         ! Calculate all basis-function hybdat%products to obtain
                         ! the overlaps with the hybdat%product-basis functions (hybdat%prodm)
                         basprod(:ng) = ( hybdat%bas1(:ng,n1,l1,itype)*hybdat%bas1(:ng,n2,l2,itype) +hybdat%bas2(:ng,n1,l1,itype)*hybdat%bas2(:ng,n2,l2,itype)) / atoms%rmsh(:ng,itype)
                         DO l = ABS(l1-l2),MIN(hybrid%lcutm1(itype),l1+l2)
                            IF(MOD(l1+l2+l,2).EQ.0) THEN
272 273
                               hybdat%nindxp1(l,itype)    = hybdat%nindxp1(l,itype) + 1
                               n                  = hybdat%nindxp1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
274 275 276 277 278 279 280
                               hybdat%prod(n,l,itype)%l1 = l1
                               hybdat%prod(n,l,itype)%l2 = l2
                               hybdat%prod(n,l,itype)%n1 = n1
                               hybdat%prod(n,l,itype)%n2 = n2
                               DO i = 1,hybrid%nindxm1(l,itype)
                                  hybdat%prodm(i,n,l,itype) = intgrf(basprod(:ng)*hybrid%basm1(:ng,i,l,itype), atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,hybdat%gridf)
                               END DO
281
                            END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
282
                         END DO
283 284

                      END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
285 286 287 288 289 290 291 292 293 294 295 296 297 298
                   END DO
                END IF
             END DO
          END DO
       END DO
       DEALLOCATE(basprod)
       CALL timestop("gen_bz and gen_wavf")


    ELSE IF (hybrid%l_hybrid ) THEN ! hybrid%l_calhf is false

       ! Reading the eig file
       !DO nk = n_start,kpts%nkpt,n_stride
       DO nk = 1,kpts%nkpt,1
Daniel Wortmann's avatar
Daniel Wortmann committed
299 300
          CALL read_eig(eig_id_hf,nk,jsp,el=el_eig, ello=ello_eig,neig=hybrid%ne_eig(nk),w_iks=results%w_iks(:,nk,jsp))
          hybrid%nobd(nk) = COUNT(results%w_iks(:hybrid%ne_eig(nk),nk,jsp) > 0.0 )
Daniel Wortmann's avatar
Daniel Wortmann committed
301
       END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
302
    
Daniel Wortmann's avatar
Daniel Wortmann committed
303
       hybrid%maxlmindx = MAXVAL((/ ( SUM( (/ (hybrid%nindx(l,itype)*(2*l+1), l=0,atoms%lmax(itype)) /) ),itype=1,atoms%ntype) /) )
Daniel Wortmann's avatar
Daniel Wortmann committed
304
       hybrid%nbands    = MIN( hybrid%bands1, DIMENSION%neigd )
Daniel Wortmann's avatar
Daniel Wortmann committed
305 306

    ENDIF ! hybrid%l_calhf
Daniel Wortmann's avatar
Daniel Wortmann committed
307 308
  END SUBROUTINE hf_setup
END MODULE m_hf_setup