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

11
   SUBROUTINE hf_setup(mpdata, hybinp, input, sym, kpts,  atoms, mpi, noco,nococonv, 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
Matthias Redies's avatar
Matthias Redies committed
18
      USE m_hybrid_core
Matthias Redies's avatar
Matthias Redies committed
19
      USE m_gen_wavf
Daniel Wortmann's avatar
Daniel Wortmann committed
20
      use m_types_hybdat
Matthias Redies's avatar
Matthias Redies committed
21 22 23

      IMPLICIT NONE

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

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

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

      ! local arrays

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

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

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

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

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

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

Matthias Redies's avatar
Matthias Redies committed
76
         if(allocated(hybdat%kveclo_eig)) deallocate(hybdat%kveclo_eig)
77
         allocate(hybdat%kveclo_eig(atoms%nlotot, kpts%nkpt), stat=ok)
78
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%kveclo_eig')
Matthias Redies's avatar
Matthias Redies committed
79
         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=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)
85 86
         END IF

Matthias Redies's avatar
Matthias Redies committed
87 88 89
         ! Reading the eig file
         DO nk = 1, kpts%nkpt
            nrec1 = kpts%nkpt*(jsp - 1) + nk
90
            CALL lapw%init(input, noco, nococonv,kpts, atoms, sym, nk, cell, sym%zrfs)
Matthias Redies's avatar
Matthias Redies committed
91
            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(input%neig*2,input%neig,noco%l_soc))
Matthias Redies's avatar
Matthias Redies committed
93
            CALL read_eig(hybdat%eig_id, 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,:input%neig) = zDebug_r(:nbasfcn,:input%neig)
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,:input%neig) = zDebug_c(:nbasfcn,:input%neig)
104 105
            END IF

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

         IF(l_exist) CLOSE(993)

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

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

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

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

155 156 157
            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
158 159 160
                  EXIT
               END IF
            END DO
161

162 163
            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
164
            END DO
165

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

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

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

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

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

         ! set up pointer pntgpt

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

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

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

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

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

         !DO nk = n_start,kpts%nkpt,n_stride
         DO nk = 1, kpts%nkpt, 1
272 273
            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
274
         END DO
275

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

279
      ENDIF ! hybdat%l_calhf
280

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

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