hf_setup.F90 13.1 KB
Newer Older
1 2 3 4 5
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------
6

Daniel Wortmann's avatar
Daniel Wortmann committed
7
MODULE m_hf_setup
8

Daniel Wortmann's avatar
Daniel Wortmann committed
9
CONTAINS
10

Matthias Redies's avatar
Matthias Redies committed
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
   SUBROUTINE hf_setup(hybrid, input, sym, kpts, DIMENSION, atoms, mpi, noco, cell, oneD, results, jsp, enpara, eig_id_hf, &
                       hybdat, it, l_real, vr0, eig_irr)
      USE m_types
      USE m_eig66_io
      USE m_util
      USE m_checkolap
      USE m_read_core
      USE m_gen_wavf

      IMPLICIT NONE

      TYPE(t_hybrid), INTENT(INOUT) :: hybrid
      TYPE(t_kpts), INTENT(IN)    :: kpts
      TYPE(t_dimension), INTENT(IN)    :: dimension
      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
      TYPE(t_enpara), INTENT(IN)    :: enpara
      TYPE(t_results), INTENT(INOUT) :: results
      TYPE(t_hybdat), INTENT(INOUT) :: hybdat

      INTEGER, INTENT(IN)    :: it
      INTEGER, INTENT(IN)    :: jsp, eig_id_hf
      REAL, INTENT(IN)    :: vr0(:, :, :)
      LOGICAL, INTENT(IN)    :: l_real

      REAL, ALLOCATABLE, INTENT(OUT)   :: eig_irr(:, :)

      ! local type variables
      TYPE(t_lapw)             :: lapw
      TYPE(t_mat), ALLOCATABLE :: zmat(:)

      ! local scalars
      INTEGER :: ok, nk, nrec1, i, j, ll, l1, l2, ng, itype, n, l, n1, n2, nn
      INTEGER :: nbasfcn
50
      LOGICAL :: l_exist
Matthias Redies's avatar
Matthias Redies committed
51 52 53 54

      ! local arrays

      REAL, ALLOCATABLE :: basprod(:)
55
      INTEGER              :: degenerat(merge(dimension%neigd*2,dimension%neigd,noco%l_soc) + 1, kpts%nkpt)
Matthias Redies's avatar
Matthias Redies committed
56 57
      LOGICAL              :: skip_kpt(kpts%nkpt)
      INTEGER              :: g(3)
58

59 60
      REAL :: zDebug_r(DIMENSION%nbasfcn,Dimension%Neigd)
      COMPLEX :: zDebug_c(DIMENSION%nbasfcn,Dimension%Neigd)
61

Matthias Redies's avatar
Matthias Redies committed
62
      skip_kpt = .FALSE.
Daniel Wortmann's avatar
Daniel Wortmann committed
63

Matthias Redies's avatar
Matthias Redies committed
64 65 66
      IF (hybrid%l_calhf) THEN
         ! Preparations for HF and hybrid functional calculation
         CALL timestart("gen_bz and gen_wavf")
67

68
         WRITE(9333,*) DIMENSION%nbasfcn, Dimension%Neigd, atoms%nlotot, kpts%nkpt
69 70
         WRITE(9333,*) ALLOCATED(hybdat%kveclo_eig)

Matthias Redies's avatar
Matthias Redies committed
71 72
         ALLOCATE (zmat(kpts%nkptf), stat=ok)
         IF (ok /= 0) STOP 'eigen_hf: failure allocation z_c'
73
         ALLOCATE (eig_irr(merge(dimension%neigd*2,dimension%neigd,noco%l_soc), kpts%nkpt), stat=ok)
Matthias Redies's avatar
Matthias Redies committed
74
         IF (ok /= 0) STOP 'eigen_hf: failure allocation eig_irr'
75
         IF(ALLOCATED(hybdat%kveclo_eig)) DEALLOCATE (hybdat%kveclo_eig) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
76 77 78 79
         ALLOCATE (hybdat%kveclo_eig(atoms%nlotot, kpts%nkpt), stat=ok)
         IF (ok /= 0) STOP 'eigen_hf: failure allocation hybdat%kveclo_eig'
         eig_irr = 0
         hybdat%kveclo_eig = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
80

81 82
         INQUIRE(file ="z",exist= l_exist)
         IF(l_exist) THEN
83 84
            IF (l_real) OPEN(unit=993,file='z',form='unformatted',access='direct',recl=DIMENSION%nbasfcn*Dimension%Neigd*8)
            IF (.NOT.l_real) OPEN(unit=993,file='z',form='unformatted',access='direct',recl=DIMENSION%nbasfcn*Dimension%Neigd*16)
85 86
         END IF

Matthias Redies's avatar
Matthias Redies committed
87 88 89 90 91
         ! Reading the eig file
         DO nk = 1, kpts%nkpt
            nrec1 = kpts%nkpt*(jsp - 1) + nk
            CALL lapw%init(input, noco, kpts, atoms, sym, nk, cell, sym%zrfs)
            nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
92
            CALL zMat(nk)%init(l_real, nbasfcn, merge(dimension%neigd*2,dimension%neigd,noco%l_soc))
Matthias Redies's avatar
Matthias Redies committed
93
            CALL read_eig(eig_id_hf, nk, jsp, zmat=zMat(nk))
94

95 96 97
            IF(l_exist.AND.zmat(1)%l_real) THEN
               READ(993,rec=nk) zDebug_r(:,:)
               zMat(nk)%data_r = 0.0
98
               zMat(nk)%data_r(:nbasfcn,:Dimension%Neigd) = zDebug_r(:nbasfcn,:Dimension%Neigd)
99 100 101 102
            END IF
            IF(l_exist.AND..NOT.zmat(1)%l_real) THEN
               READ(993,rec=nk) zDebug_c(:,:)
               zMat(nk)%data_c = 0.0
103
               zMat(nk)%data_c(:nbasfcn,:Dimension%Neigd) = zDebug_c(:nbasfcn,:Dimension%Neigd)
104 105 106 107 108
            END IF

            WRITE(9333,*) SHAPE(eig_irr)
            WRITE(9333,*) SHAPE(results%eig)

Matthias Redies's avatar
Matthias Redies committed
109 110 111
            eig_irr(:, nk) = results%eig(:, nk, jsp)
            hybrid%ne_eig(nk) = results%neig(nk, jsp)
         END DO
112 113 114

         IF(l_exist) CLOSE(993)

Matthias Redies's avatar
Matthias Redies committed
115 116 117
         !Allocate further space
         DO nk = kpts%nkpt + 1, kpts%nkptf
            nbasfcn = zMat(kpts%bkp(nk))%matsize1
118
            CALL zMat(nk)%init(l_real, nbasfcn, merge(dimension%neigd*2,dimension%neigd,noco%l_soc))
Matthias Redies's avatar
Matthias Redies committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
         END DO

         !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
         hybrid%nobd = 0
         DO nk = 1, kpts%nkpt
            DO i = 1, hybrid%ne_eig(nk)
               DO j = i + 1, hybrid%ne_eig(nk)
                  IF (ABS(results%eig(i, nk, jsp) - results%eig(j, nk, jsp)) < 1E-07) THEN !0.015
                     degenerat(i, nk) = degenerat(i, nk) + 1
                  END IF
               END DO
140 141
            END DO

Matthias Redies's avatar
Matthias Redies committed
142 143 144
            DO i = 1, hybrid%ne_eig(nk)
               IF ((degenerat(i, nk) /= 1) .OR. (degenerat(i, nk) /= 0)) degenerat(i + 1:i + degenerat(i, nk) - 1, nk) = 0
            END DO
145

Matthias Redies's avatar
Matthias Redies committed
146
            ! set the size of the exchange matrix in the space of the wavefunctions
147

Matthias Redies's avatar
Matthias Redies committed
148 149 150 151 152 153 154 155
            hybrid%nbands(nk) = hybrid%bands1
            IF (hybrid%nbands(nk) > hybrid%ne_eig(nk)) THEN
               IF (mpi%irank == 0) THEN
                  WRITE (*, *) ' maximum for hybrid%nbands is', hybrid%ne_eig(nk)
                  WRITE (*, *) ' increase energy window to obtain enough eigenvalues'
                  WRITE (*, *) ' set hybrid%nbands equal to hybrid%ne_eig'
               END IF
               hybrid%nbands(nk) = hybrid%ne_eig(nk)
156 157
            END IF

Matthias Redies's avatar
Matthias Redies committed
158 159 160 161 162 163
            DO i = hybrid%nbands(nk) - 1, 1, -1
               IF ((degenerat(i, nk) >= 1) .AND. (degenerat(i, nk) + i - 1 /= hybrid%nbands(nk))) THEN
                  hybrid%nbands(nk) = i + degenerat(i, nk) - 1
                  EXIT
               END IF
            END DO
164

Matthias Redies's avatar
Matthias Redies committed
165 166 167
            DO i = 1, hybrid%ne_eig(nk)
               IF (results%w_iks(i, nk, jsp) > 0.0) hybrid%nobd(nk) = hybrid%nobd(nk) + 1
            END DO
168

Matthias Redies's avatar
Matthias Redies committed
169 170 171 172 173 174 175 176 177
            IF (hybrid%nobd(nk) > hybrid%nbands(nk)) THEN
               WRITE (*, *) 'k-point: ', nk
               WRITE (*, *) 'number of bands:          ', hybrid%nbands(nk)
               WRITE (*, *) 'number of occupied bands: ', hybrid%nobd(nk)
               CALL judft_warn("More occupied bands than total no of bands!?")
               hybrid%nbands(nk) = hybrid%nobd(nk)
            END IF
            PRINT *, "bands:", nk, hybrid%nobd(nk), hybrid%nbands(nk), hybrid%ne_eig(nk)
         END DO
178

Matthias Redies's avatar
Matthias Redies committed
179 180 181 182 183
         ! spread hybrid%nobd from IBZ to whole BZ
         DO nk = 1, kpts%nkptf
            i = kpts%bkp(nk)
            hybrid%nobd(nk) = hybrid%nobd(i)
         END DO
184

Matthias Redies's avatar
Matthias Redies committed
185 186 187
         ! generate eigenvectors z and MT coefficients from the previous iteration at all k-points
         CALL gen_wavf(kpts%nkpt, kpts, it, sym, atoms, enpara%el0(:, :, jsp), enpara%ello0(:, :, jsp), cell, dimension, &
                       hybrid, vr0, hybdat, noco, oneD, mpi, input, jsp, zmat)
188

Matthias Redies's avatar
Matthias Redies committed
189 190 191
         ! 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)
192

Matthias Redies's avatar
Matthias Redies committed
193 194 195 196 197 198 199
         ! check olap between core-basis/core-valence/basis-basis
         CALL checkolap(atoms, hybdat, hybrid, kpts%nkpt, kpts, dimension, mpi, skip_kpt, &
                        input, sym, noco, cell, lapw, jsp)

         ! set up pointer pntgpt

         ! setup dimension of pntgpt
200
         IF(ALLOCATED(hybdat%pntgptd)) DEALLOCATE(hybdat%pntgptd) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
201 202 203 204 205 206 207 208 209
         ALLOCATE (hybdat%pntgptd(3))
         hybdat%pntgptd = 0
         DO nk = 1, kpts%nkptf
            CALL lapw%init(input, noco, kpts, atoms, sym, nk, cell, sym%zrfs)
            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

210
         IF(ALLOCATED(hybdat%pntgpt)) DEALLOCATE(hybdat%pntgpt) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
211 212 213 214 215 216 217 218 219 220
         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 /= 0) STOP 'eigen_hf: failure allocation pntgpt'
         hybdat%pntgpt = 0
         DO nk = 1, kpts%nkptf
            CALL lapw%init(input, noco, kpts, atoms, sym, nk, cell, sym%zrfs)
            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
221
         END DO
Matthias Redies's avatar
Matthias Redies committed
222 223 224

         ALLOCATE (basprod(atoms%jmtd), stat=ok)
         IF (ok /= 0) STOP 'eigen_hf: failure allocation basprod'
225
         IF(ALLOCATED(hybdat%prodm)) DEALLOCATE(hybdat%prodm) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
226 227
         ALLOCATE (hybdat%prodm(hybrid%maxindxm1, hybrid%maxindxp1, 0:hybrid%maxlcutm1, atoms%ntype), stat=ok)
         IF (ok /= 0) STOP 'eigen_hf: failure allocation hybdat%prodm'
228
         IF(ALLOCATED(hybdat%prod)) DEALLOCATE(hybdat%prod) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
229 230 231 232
         ALLOCATE (hybdat%prod(hybrid%maxindxp1, 0:hybrid%maxlcutm1, atoms%ntype), stat=ok)
         IF (ok /= 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
233
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE(hybdat%nindxp1) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
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
         ALLOCATE (hybdat%nindxp1(0:hybrid%maxlcutm1, atoms%ntype))
         hybdat%nindxp1 = 0
         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) <= hybrid%lcutm1(itype)) THEN
                     DO n2 = 1, hybrid%nindx(l2, itype)
                        nn = hybrid%nindx(l1, itype)
                        IF (l1 == 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) == 0) THEN
                                 hybdat%nindxp1(l, itype) = hybdat%nindxp1(l, itype) + 1
                                 n = hybdat%nindxp1(l, itype)
                                 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
                              END IF
                           END DO
264 265
                        END DO
                     END DO
Matthias Redies's avatar
Matthias Redies committed
266 267
                  END IF
               END DO
268 269
            END DO
         END DO
Matthias Redies's avatar
Matthias Redies committed
270 271
         DEALLOCATE (basprod)
         CALL timestop("gen_bz and gen_wavf")
272

Matthias Redies's avatar
Matthias Redies committed
273 274 275 276 277 278 279
      ELSE IF (hybrid%l_hybrid) THEN ! hybrid%l_calhf is false

         !DO nk = n_start,kpts%nkpt,n_stride
         DO nk = 1, kpts%nkpt, 1
            hybrid%ne_eig(nk) = results%neig(nk, jsp)
            hybrid%nobd(nk) = COUNT(results%w_iks(:hybrid%ne_eig(nk), nk, jsp) > 0.0)
         END DO
280

Matthias Redies's avatar
Matthias Redies committed
281 282
         hybrid%maxlmindx = MAXVAL((/(SUM((/(hybrid%nindx(l, itype)*(2*l + 1), l=0, atoms%lmax(itype))/)), itype=1, atoms%ntype)/))
         hybrid%nbands = MIN(hybrid%bands1, DIMENSION%neigd)
283

Matthias Redies's avatar
Matthias Redies committed
284
      ENDIF ! hybrid%l_calhf
285

Matthias Redies's avatar
Matthias Redies committed
286
   END SUBROUTINE hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
287

Daniel Wortmann's avatar
Daniel Wortmann committed
288
END MODULE m_hf_setup