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 12
   SUBROUTINE hf_setup(mpdata, hybinp, input, sym, kpts,  atoms, mpi, noco,nococonv, &
                       cell, oneD, results, jsp, enpara, &
13
                       hybdat, l_real, vr0, eig_irr)
Matthias Redies's avatar
Matthias Redies committed
14 15 16
      USE m_types
      USE m_eig66_io
      USE m_util
Matthias Redies's avatar
Matthias Redies committed
17
      USE m_intgrf
Matthias Redies's avatar
Matthias Redies committed
18
      USE m_checkolap
Matthias Redies's avatar
Matthias Redies committed
19
      USE m_hybrid_core
Matthias Redies's avatar
Matthias Redies committed
20
      USE m_gen_wavf
Daniel Wortmann's avatar
Daniel Wortmann committed
21
      use m_types_hybdat
Matthias Redies's avatar
Matthias Redies committed
22 23 24

      IMPLICIT NONE

25
      TYPE(t_mpdata), INTENT(inout)   :: mpdata
26
      TYPE(t_hybinp), INTENT(IN) :: hybinp
Matthias Redies's avatar
Matthias Redies committed
27 28 29 30
      TYPE(t_kpts), INTENT(IN)    :: kpts
      TYPE(t_atoms), INTENT(IN)    :: atoms
      TYPE(t_mpi), INTENT(IN)    :: mpi
      TYPE(t_noco), INTENT(IN)    :: noco
31
      TYPE(t_nococonv), INTENT(IN)    :: nococonv
Matthias Redies's avatar
Matthias Redies committed
32 33 34 35 36 37 38 39
      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

Matthias Redies's avatar
Matthias Redies committed
40
      INTEGER, INTENT(IN)    :: jsp
Matthias Redies's avatar
Matthias Redies committed
41 42 43
      REAL, INTENT(IN)    :: vr0(:, :, :)
      LOGICAL, INTENT(IN)    :: l_real

Matthias Redies's avatar
Matthias Redies committed
44
      REAL, ALLOCATABLE, INTENT(INOUT)   :: eig_irr(:, :)
Matthias Redies's avatar
Matthias Redies committed
45 46 47 48 49 50 51

      ! 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
52
      INTEGER :: nbasfcn, n_dim
53
      LOGICAL :: l_exist
Matthias Redies's avatar
Matthias Redies committed
54 55 56 57

      ! local arrays

      REAL, ALLOCATABLE :: basprod(:)
58
      INTEGER              :: degenerat(merge(input%neig*2,input%neig,noco%l_soc) + 1, kpts%nkpt)
Matthias Redies's avatar
Matthias Redies committed
59
      LOGICAL              :: skip_kpt(kpts%nkpt)
60

61 62
      REAL :: zDebug_r(lapw_dim_nbasfcn,input%neig)
      COMPLEX :: zDebug_c(lapw_dim_nbasfcn,input%neig)
63

Matthias Redies's avatar
Matthias Redies committed
64
      skip_kpt = .FALSE.
Daniel Wortmann's avatar
Daniel Wortmann committed
65

66
      IF (hybdat%l_calhf) THEN
Matthias Redies's avatar
Matthias Redies committed
67
         ! Preparations for HF and hybinp functional calculation
Matthias Redies's avatar
Matthias Redies committed
68
         CALL timestart("gen_bz and gen_wavf")
69

70
         allocate(zmat(kpts%nkptf), stat=ok)
71
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation z_c')
72 73

         if(.not. allocated(eig_irr)) allocate(eig_irr(input%neig, kpts%nkpt), stat=ok)
74
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation eig_irr')
75
         eig_irr = 0
Matthias Redies's avatar
Matthias Redies committed
76

Matthias Redies's avatar
Matthias Redies committed
77
         if(allocated(hybdat%kveclo_eig)) deallocate(hybdat%kveclo_eig)
78
         allocate(hybdat%kveclo_eig(atoms%nlotot, kpts%nkpt), stat=ok)
79
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%kveclo_eig')
Matthias Redies's avatar
Matthias Redies committed
80
         hybdat%kveclo_eig = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
81

82 83
         INQUIRE(file ="z",exist= l_exist)
         IF(l_exist) THEN
84 85
            IF (l_real) OPEN(unit=993,file='z',form='unformatted',access='direct',recl=lapw%dim_nbasfcn()*input%neig*8)
            IF (.NOT.l_real) OPEN(unit=993,file='z',form='unformatted',access='direct',recl=lapw%dim_nbasfcn()*input%neig*16)
86 87
         END IF

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

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

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

         IF(l_exist) CLOSE(993)

Matthias Redies's avatar
Matthias Redies committed
113 114 115
         !Allocate further space
         DO nk = kpts%nkpt + 1, kpts%nkptf
            nbasfcn = zMat(kpts%bkp(nk))%matsize1
116
            CALL zMat(nk)%init(l_real, nbasfcn, merge(input%neig*2,input%neig,noco%l_soc))
Matthias Redies's avatar
Matthias Redies committed
117 118 119 120 121 122 123 124 125 126 127 128 129
         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
130
         hybdat%nobd(:,jsp) = 0
Matthias Redies's avatar
Matthias Redies committed
131
         DO nk = 1, kpts%nkpt
132 133
            DO i = 1, hybdat%ne_eig(nk)
               DO j = i + 1, hybdat%ne_eig(nk)
Matthias Redies's avatar
Matthias Redies committed
134 135 136 137
                  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
138 139
            END DO

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

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

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

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

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

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

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

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

Matthias Redies's avatar
Matthias Redies committed
188
         ! generate core wave functions (-> core1/2(jmtd,hybdat%nindxc,0:lmaxc,ntype) )
189
         CALL corewf(atoms, jsp, input,  vr0, hybdat%lmaxcd, hybdat%maxindxc, mpi, &
Matthias Redies's avatar
Matthias Redies committed
190
                     hybdat%lmaxc, hybdat%nindxc, hybdat%core1, hybdat%core2, hybdat%eig_c)
191

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

         ! set up pointer pntgpt

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

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

221
         allocate(basprod(atoms%jmtd), stat=ok)
222
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation basprod')
Matthias Redies's avatar
Matthias Redies committed
223
         IF(ALLOCATED(hybdat%prodm)) DEALLOCATE(hybdat%prodm)
224
         allocate(hybdat%prodm(maxval(mpdata%num_radbasfn), mpdata%max_indx_p_1, 0:maxval(hybinp%lcutm1), atoms%ntype), stat=ok)
225
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%prodm')
226

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

269
      ELSE IF (hybinp%l_hybrid) THEN ! hybdat%l_calhf is false
Matthias Redies's avatar
Matthias Redies committed
270 271 272

         !DO nk = n_start,kpts%nkpt,n_stride
         DO nk = 1, kpts%nkpt, 1
273 274
            hybdat%ne_eig(nk) = results%neig(nk, jsp)
            hybdat%nobd(nk,jsp) = COUNT(results%w_iks(:hybdat%ne_eig(nk), nk, jsp) > 0.0)
Matthias Redies's avatar
Matthias Redies committed
275
         END DO
276

277
         hybdat%maxlmindx = MAXVAL([(SUM([(mpdata%num_radfun_per_l(l, itype)*(2*l + 1), l=0, atoms%lmax(itype))]), itype=1, atoms%ntype)])
278
         hybdat%nbands = MIN(hybinp%bands1, input%neig)
279

280
      ENDIF ! hybdat%l_calhf
281

Matthias Redies's avatar
Matthias Redies committed
282
   END SUBROUTINE hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
283

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