hf_setup.F90 11 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
   SUBROUTINE hf_setup(hybrid, input, sym, kpts, DIMENSION, atoms, mpi, noco, cell, oneD, results, jsp, enpara, eig_id_hf, &
12
                       hybdat, l_real, vr0, eig_irr)
Matthias Redies's avatar
Matthias Redies committed
13 14 15 16
      USE m_types
      USE m_eig66_io
      USE m_util
      USE m_checkolap
17
      USE m_hybrid_core
Matthias Redies's avatar
Matthias Redies committed
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
      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)    :: 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
Matthias Redies's avatar
Matthias Redies committed
48
      INTEGER :: nbasfcn, n_dim
Matthias Redies's avatar
Matthias Redies committed
49 50 51 52 53 54

      ! local arrays

      REAL, ALLOCATABLE :: basprod(:)
      INTEGER              :: degenerat(DIMENSION%neigd2 + 1, kpts%nkpt)
      LOGICAL              :: skip_kpt(kpts%nkpt)
55

Matthias Redies's avatar
Matthias Redies committed
56
      skip_kpt = .FALSE.
Daniel Wortmann's avatar
Daniel Wortmann committed
57

Matthias Redies's avatar
Matthias Redies committed
58 59 60
      IF (hybrid%l_calhf) THEN
         ! Preparations for HF and hybrid functional calculation
         CALL timestart("gen_bz and gen_wavf")
61

Matthias Redies's avatar
Matthias Redies committed
62
         ALLOCATE (zmat(kpts%nkptf), stat=ok)
63
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation z_c')
Matthias Redies's avatar
Matthias Redies committed
64
         ALLOCATE (eig_irr(DIMENSION%neigd2, kpts%nkpt), stat=ok)
65
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation eig_irr')
Matthias Redies's avatar
Matthias Redies committed
66
         ALLOCATE (hybdat%kveclo_eig(atoms%nlotot, kpts%nkpt), stat=ok)
67
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%kveclo_eig')
Matthias Redies's avatar
Matthias Redies committed
68 69
         eig_irr = 0
         hybdat%kveclo_eig = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
70

Matthias Redies's avatar
Matthias Redies committed
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
         ! 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)
            CALL zMat(nk)%init(l_real, nbasfcn, dimension%neigd2)
            CALL read_eig(eig_id_hf, nk, jsp, zmat=zMat(nk))
            eig_irr(:, nk) = results%eig(:, nk, jsp)
            hybrid%ne_eig(nk) = results%neig(nk, jsp)
         END DO
         !Allocate further space
         DO nk = kpts%nkpt + 1, kpts%nkptf
            nbasfcn = zMat(kpts%bkp(nk))%matsize1
            CALL zMat(nk)%init(l_real, nbasfcn, dimension%neigd2)
         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
106 107
            END DO

Matthias Redies's avatar
Matthias Redies committed
108 109 110
            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
111

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

Matthias Redies's avatar
Matthias Redies committed
114 115 116 117 118 119 120 121
            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)
122 123
            END IF

Matthias Redies's avatar
Matthias Redies committed
124 125 126 127 128 129
            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
130

Matthias Redies's avatar
Matthias Redies committed
131 132 133
            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
134

Matthias Redies's avatar
Matthias Redies committed
135 136 137 138 139 140 141 142 143
            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
144

Matthias Redies's avatar
Matthias Redies committed
145 146 147 148 149
         ! 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
150

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

Matthias Redies's avatar
Matthias Redies committed
155 156 157
         ! 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)
158

Matthias Redies's avatar
Matthias Redies committed
159
         ! check olap between core-basis/core-valence/basis-basis
Matthias Redies's avatar
Matthias Redies committed
160
         CALL checkolap(atoms, hybdat, hybrid, kpts%nkpt, kpts, dimension, mpi, &
Matthias Redies's avatar
Matthias Redies committed
161 162 163 164 165 166 167 168 169
                        input, sym, noco, cell, lapw, jsp)

         ! set up pointer pntgpt

         ! setup dimension of pntgpt
         ALLOCATE (hybdat%pntgptd(3))
         hybdat%pntgptd = 0
         DO nk = 1, kpts%nkptf
            CALL lapw%init(input, noco, kpts, atoms, sym, nk, cell, sym%zrfs)
Matthias Redies's avatar
Matthias Redies committed
170 171 172
            do n_dim = 1,3
               hybdat%pntgptd(n_dim) = MAXVAL([(ABS(lapw%gvec(n_dim,i, jsp)), i=1), lapw%nv(jsp), hybdat%pntgptd(n_dim) ])
            end do
Matthias Redies's avatar
Matthias Redies committed
173 174 175 176
         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)
177
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation pntgpt')
Matthias Redies's avatar
Matthias Redies committed
178 179 180 181
         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)
182
               hybdat%pntgpt(lapw%gvec(1,i,jsp), lapw%gvec(2,i,jsp), lapw%gvec(3,i,jsp), nk) = i
Matthias Redies's avatar
Matthias Redies committed
183
            END DO
184
         END DO
Matthias Redies's avatar
Matthias Redies committed
185 186

         ALLOCATE (basprod(atoms%jmtd), stat=ok)
187
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation basprod')
Matthias Redies's avatar
Matthias Redies committed
188
         ALLOCATE (hybdat%prodm(hybrid%maxindxm1, hybrid%maxindxp1, 0:hybrid%maxlcutm1, atoms%ntype), stat=ok)
189
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%prodm')
Matthias Redies's avatar
Matthias Redies committed
190
         ALLOCATE (hybdat%prod(hybrid%maxindxp1, 0:hybrid%maxlcutm1, atoms%ntype), stat=ok)
191
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%prod')
Matthias Redies's avatar
Matthias Redies committed
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
         basprod = 0; hybdat%prodm = 0; hybdat%prod%l1 = 0; hybdat%prod%l2 = 0
         hybdat%prod%n1 = 0; hybdat%prod%n2 = 0
         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
224 225
                        END DO
                     END DO
Matthias Redies's avatar
Matthias Redies committed
226 227
                  END IF
               END DO
228 229
            END DO
         END DO
Matthias Redies's avatar
Matthias Redies committed
230 231
         DEALLOCATE (basprod)
         CALL timestop("gen_bz and gen_wavf")
232

Matthias Redies's avatar
Matthias Redies committed
233 234 235 236 237 238 239
      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
240

Matthias Redies's avatar
Matthias Redies committed
241 242
         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)
243

Matthias Redies's avatar
Matthias Redies committed
244
      ENDIF ! hybrid%l_calhf
245

Matthias Redies's avatar
Matthias Redies committed
246
   END SUBROUTINE hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
247

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