hf_setup.F90 12.5 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(mpbasis, 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
      USE m_types
      USE m_eig66_io
      USE m_util
Matthias Redies's avatar
Matthias Redies committed
16
      USE m_intgrf
Matthias Redies's avatar
Matthias Redies committed
17
      USE m_checkolap
18
      USE m_hybrid_core
Matthias Redies's avatar
Matthias Redies committed
19 20 21 22
      USE m_gen_wavf

      IMPLICIT NONE

Matthias Redies's avatar
Matthias Redies committed
23
      TYPE(t_mpbasis), INTENT(in)   :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
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
      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
50
      INTEGER :: nbasfcn, n_dim
51
      LOGICAL :: l_exist
Matthias Redies's avatar
Matthias Redies committed
52 53 54 55 56 57

      ! local arrays

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

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

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
         allocate(zmat(kpts%nkptf), stat=ok)
69
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation z_c')
70
         allocate(eig_irr(DIMENSION%neigd2, kpts%nkpt), stat=ok)
71
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation eig_irr')
72
         allocate(hybdat%kveclo_eig(atoms%nlotot, kpts%nkpt), stat=ok)
73
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%kveclo_eig')
Matthias Redies's avatar
Matthias Redies committed
74 75
         eig_irr = 0
         hybdat%kveclo_eig = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
76

77 78 79 80 81 82
         INQUIRE(file ="z",exist= l_exist)
         IF(l_exist) THEN
            IF (l_real) OPEN(unit=993,file='z',form='unformatted',access='direct',recl=DIMENSION%nbasfcn*DIMENSION%neigd2*8)
            IF (.NOT.l_real) OPEN(unit=993,file='z',form='unformatted',access='direct',recl=DIMENSION%nbasfcn*DIMENSION%neigd2*16)
         END IF

Matthias Redies's avatar
Matthias Redies committed
83 84 85 86 87 88 89
         ! 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))
Matthias Redies's avatar
Matthias Redies committed
90

91 92 93 94 95 96 97 98 99 100 101 102 103 104
            IF(l_exist.AND.zmat(1)%l_real) THEN
               READ(993,rec=nk) zDebug_r(:,:)
               zMat(nk)%data_r = 0.0
               zMat(nk)%data_r(:nbasfcn,:DIMENSION%neigd2) = zDebug_r(:nbasfcn,:DIMENSION%neigd2)
            END IF
            IF(l_exist.AND..NOT.zmat(1)%l_real) THEN
               READ(993,rec=nk) zDebug_c(:,:)
               zMat(nk)%data_c = 0.0
               zMat(nk)%data_c(:nbasfcn,:DIMENSION%neigd2) = zDebug_c(:nbasfcn,:DIMENSION%neigd2)
            END IF

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

Matthias Redies's avatar
Matthias Redies committed
105 106 107
            eig_irr(:, nk) = results%eig(:, nk, jsp)
            hybrid%ne_eig(nk) = results%neig(nk, jsp)
         END DO
108 109 110

         IF(l_exist) CLOSE(993)

Matthias Redies's avatar
Matthias Redies committed
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
         !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
128
         hybrid%nobd(:,jsp) = 0
Matthias Redies's avatar
Matthias Redies committed
129 130 131 132 133 134 135
         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
136 137
            END DO

Matthias Redies's avatar
Matthias Redies committed
138 139 140
            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
141

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

Matthias Redies's avatar
Matthias Redies committed
144 145 146 147 148 149 150 151
            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)
152 153
            END IF

Matthias Redies's avatar
Matthias Redies committed
154 155 156 157 158 159
            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
160

Matthias Redies's avatar
Matthias Redies committed
161
            DO i = 1, hybrid%ne_eig(nk)
162
               IF (results%w_iks(i, nk, jsp) > 0.0) hybrid%nobd(nk,jsp) = hybrid%nobd(nk,jsp) + 1
Matthias Redies's avatar
Matthias Redies committed
163
            END DO
164

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

Matthias Redies's avatar
Matthias Redies committed
175 176 177
         ! spread hybrid%nobd from IBZ to whole BZ
         DO nk = 1, kpts%nkptf
            i = kpts%bkp(nk)
178
            hybrid%nobd(nk,jsp) = hybrid%nobd(i,jsp)
Matthias Redies's avatar
Matthias Redies committed
179
         END DO
180

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

Matthias Redies's avatar
Matthias Redies committed
185 186 187
         ! 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)
188

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

         ! set up pointer pntgpt

         ! setup dimension of pntgpt
196
         IF(ALLOCATED(hybdat%pntgptd)) DEALLOCATE(hybdat%pntgptd) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
197 198 199 200
         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
201
            do n_dim = 1,3
Matthias Redies's avatar
again  
Matthias Redies committed
202
               hybdat%pntgptd(n_dim) = MAXVAL([(ABS(lapw%gvec(n_dim,i,jsp)), i=1, lapw%nv(jsp)), hybdat%pntgptd(n_dim)])
Matthias Redies's avatar
Matthias Redies committed
203
            end do
Matthias Redies's avatar
Matthias Redies committed
204 205
         END DO

206
         IF(ALLOCATED(hybdat%pntgpt)) DEALLOCATE(hybdat%pntgpt) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
207 208
         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)
209
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation pntgpt')
Matthias Redies's avatar
Matthias Redies committed
210 211 212 213
         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)
214
               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
215
            END DO
216
         END DO
Matthias Redies's avatar
Matthias Redies committed
217

218
         allocate(basprod(atoms%jmtd), stat=ok)
219
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation basprod')
220
         allocate(hybdat%prodm(maxval(mpbasis%num_rad_bas_fun), hybrid%max_indx_p_1, 0:maxval(hybrid%lcutm1), atoms%ntype), stat=ok)
221
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%prodm')
222

223
         call hybdat%prod%init(hybrid, atoms)
224

Matthias Redies's avatar
Matthias Redies committed
225 226
         basprod = 0; hybdat%prodm = 0; hybdat%prod%l1 = 0; hybdat%prod%l2 = 0
         hybdat%prod%n1 = 0; hybdat%prod%n2 = 0
227
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE(hybdat%nindxp1) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
228
         ALLOCATE (hybdat%nindxp1(0:maxval(hybrid%lcutm1), atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
229 230 231 232 233 234 235
         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
236 237
                     DO n2 = 1, hybrid%num_radfun_per_l(l2, itype)
                        nn = hybrid%num_radfun_per_l(l1, itype)
Matthias Redies's avatar
Matthias Redies committed
238 239 240 241 242 243 244 245 246 247
                        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)
248 249 250 251
                                 hybdat%prod%l1(n,l,itype) = l1
                                 hybdat%prod%l2(n,l,itype) = l2
                                 hybdat%prod%n1(n,l,itype) = n1
                                 hybdat%prod%n2(n,l,itype) = n2
252
                                 DO i = 1, mpbasis%num_rad_bas_fun(l, itype)
Matthias Redies's avatar
Matthias Redies committed
253
                                    hybdat%prodm(i, n, l, itype) = intgrf(basprod(:ng)*mpbasis%radbasfn_mt(:ng, i, l, itype), atoms%jri, &
Matthias Redies's avatar
Matthias Redies committed
254 255 256 257
                                                                          atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
                                 END DO
                              END IF
                           END DO
258 259
                        END DO
                     END DO
Matthias Redies's avatar
Matthias Redies committed
260 261
                  END IF
               END DO
262 263
            END DO
         END DO
264
         deallocate(basprod)
Matthias Redies's avatar
Matthias Redies committed
265
         CALL timestop("gen_bz and gen_wavf")
266

Matthias Redies's avatar
Matthias Redies committed
267 268 269 270 271
      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)
272
            hybrid%nobd(nk,jsp) = COUNT(results%w_iks(:hybrid%ne_eig(nk), nk, jsp) > 0.0)
Matthias Redies's avatar
Matthias Redies committed
273
         END DO
274

Matthias Redies's avatar
Matthias Redies committed
275
         hybrid%maxlmindx = MAXVAL([(SUM([(hybrid%num_radfun_per_l(l, itype)*(2*l + 1), l=0, atoms%lmax(itype))]), itype=1, atoms%ntype)])
Matthias Redies's avatar
Matthias Redies committed
276
         hybrid%nbands = MIN(hybrid%bands1, DIMENSION%neigd)
277

Matthias Redies's avatar
Matthias Redies committed
278
      ENDIF ! hybrid%l_calhf
279

Matthias Redies's avatar
Matthias Redies committed
280
   END SUBROUTINE hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
281

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