hf_setup.F90 11.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
   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

      ! local type variables
      TYPE(t_lapw)             :: lapw

      ! 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)
58

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

62
      IF (hybdat%l_calhf) THEN
Matthias Redies's avatar
Matthias Redies committed
63
         ! Preparations for HF and hybinp functional calculation
Matthias Redies's avatar
Matthias Redies committed
64
         CALL timestart("gen_bz and gen_wavf")
65
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation z_c')
66 67

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

Matthias Redies's avatar
Matthias Redies committed
71
         if(allocated(hybdat%kveclo_eig)) deallocate(hybdat%kveclo_eig)
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
         hybdat%kveclo_eig = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
75

76

Matthias Redies's avatar
Matthias Redies committed
77
         ! Reading the eig file
Matthias Redies's avatar
Matthias Redies committed
78
         call timestart("eig stuff")
Matthias Redies's avatar
Matthias Redies committed
79 80
         DO nk = 1, kpts%nkpt
            nrec1 = kpts%nkpt*(jsp - 1) + nk
81
            CALL lapw%init(input, noco, nococonv,kpts, atoms, sym, nk, cell, sym%zrfs)
Matthias Redies's avatar
Matthias Redies committed
82
            nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
83

Matthias Redies's avatar
Matthias Redies committed
84
            eig_irr(:, nk) = results%eig(:, nk, jsp)
85
            hybdat%ne_eig(nk) = results%neig(nk, jsp)
Matthias Redies's avatar
Matthias Redies committed
86
         END DO
Matthias Redies's avatar
Matthias Redies committed
87
         call timestop("eig stuff")
Matthias Redies's avatar
Matthias Redies committed
88 89 90 91 92 93 94 95 96 97 98 99

         !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
100
         hybdat%nobd(:,jsp) = 0
Matthias Redies's avatar
Matthias Redies committed
101
         DO nk = 1, kpts%nkpt
102 103
            DO i = 1, hybdat%ne_eig(nk)
               DO j = i + 1, hybdat%ne_eig(nk)
Matthias Redies's avatar
Matthias Redies committed
104 105 106 107
                  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
108 109
            END DO

110
            DO i = 1, hybdat%ne_eig(nk)
Matthias Redies's avatar
Matthias Redies committed
111 112
               IF ((degenerat(i, nk) /= 1) .OR. (degenerat(i, nk) /= 0)) degenerat(i + 1:i + degenerat(i, nk) - 1, nk) = 0
            END DO
113

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

116 117
            hybdat%nbands(nk) = hybinp%bands1
            IF (hybdat%nbands(nk) > hybdat%ne_eig(nk)) THEN
Matthias Redies's avatar
Matthias Redies committed
118
               IF (mpi%irank == 0) THEN
119
                  WRITE (*, *) ' maximum for hybdat%nbands is', hybdat%ne_eig(nk)
Matthias Redies's avatar
Matthias Redies committed
120
                  WRITE (*, *) ' increase energy window to obtain enough eigenvalues'
121
                  WRITE (*, *) ' set hybdat%nbands equal to hybdat%ne_eig'
Matthias Redies's avatar
Matthias Redies committed
122
               END IF
123
               hybdat%nbands(nk) = hybdat%ne_eig(nk)
124 125
            END IF

126 127 128
            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
129 130 131
                  EXIT
               END IF
            END DO
132

133 134
            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
135
            END DO
136

137
            IF (hybdat%nobd(nk,jsp) > hybdat%nbands(nk)) THEN
Matthias Redies's avatar
Matthias Redies committed
138
               WRITE (*, *) 'k-point: ', nk
139 140
               WRITE (*, *) 'number of bands:          ', hybdat%nbands(nk)
               WRITE (*, *) 'number of occupied bands: ', hybdat%nobd(nk,jsp)
Matthias Redies's avatar
Matthias Redies committed
141
               CALL judft_warn("More occupied bands than total no of bands!?")
142
               hybdat%nbands(nk) = hybdat%nobd(nk,jsp)
Matthias Redies's avatar
Matthias Redies committed
143 144
            END IF
         END DO
145

146
         ! spread hybdat%nobd from IBZ to whole BZ
Matthias Redies's avatar
Matthias Redies committed
147 148
         DO nk = 1, kpts%nkptf
            i = kpts%bkp(nk)
Matthias Redies's avatar
Matthias Redies committed
149
            hybdat%nbands(nk) = hybdat%nbands(i)
150
            hybdat%nobd(nk,jsp) = hybdat%nobd(i,jsp)
Matthias Redies's avatar
Matthias Redies committed
151
         END DO
152

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

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

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

         ! set up pointer pntgpt

         ! setup dimension of pntgpt
168
         IF(ALLOCATED(hybdat%pntgptd)) DEALLOCATE(hybdat%pntgptd) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
169 170 171
         ALLOCATE (hybdat%pntgptd(3))
         hybdat%pntgptd = 0
         DO nk = 1, kpts%nkptf
172
            CALL lapw%init(input, noco, nococonv,kpts, atoms, sym, nk, cell, sym%zrfs)
Matthias Redies's avatar
Matthias Redies committed
173
            do n_dim = 1,3
Matthias Redies's avatar
again  
Matthias Redies committed
174
               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
175
            end do
Matthias Redies's avatar
Matthias Redies committed
176 177
         END DO

178
         IF(ALLOCATED(hybdat%pntgpt)) DEALLOCATE(hybdat%pntgpt) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
179 180
         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)
181
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation pntgpt')
Matthias Redies's avatar
Matthias Redies committed
182 183
         hybdat%pntgpt = 0
         DO nk = 1, kpts%nkptf
184
            CALL lapw%init(input, noco, nococonv,kpts, atoms, sym, nk, cell, sym%zrfs)
Matthias Redies's avatar
Matthias Redies committed
185
            DO i = 1, lapw%nv(jsp)
186
               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
187
            END DO
188
         END DO
Matthias Redies's avatar
Matthias Redies committed
189

190
         allocate(basprod(atoms%jmtd), stat=ok)
191
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation basprod')
Matthias Redies's avatar
Matthias Redies committed
192
         IF(ALLOCATED(hybdat%prodm)) DEALLOCATE(hybdat%prodm)
193
         allocate(hybdat%prodm(maxval(mpdata%num_radbasfn), mpdata%max_indx_p_1, 0:maxval(hybinp%lcutm1), atoms%ntype), stat=ok)
194
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%prodm')
195

196 197
         basprod = 0; hybdat%prodm = 0; mpdata%l1 = 0; mpdata%l2 = 0
         mpdata%n1 = 0; mpdata%n2 = 0
198
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE(hybdat%nindxp1) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
199
         ALLOCATE (hybdat%nindxp1(0:maxval(hybinp%lcutm1), atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
200 201 202
         hybdat%nindxp1 = 0
         DO itype = 1, atoms%ntype
            ng = atoms%jri(itype)
Matthias Redies's avatar
Matthias Redies committed
203
            DO l2 = 0, MIN(atoms%lmax(itype), hybinp%lcutwf(itype))
Matthias Redies's avatar
Matthias Redies committed
204 205
               ll = l2
               DO l1 = 0, ll
Matthias Redies's avatar
Matthias Redies committed
206
                  IF (ABS(l1 - l2) <= hybinp%lcutm1(itype)) THEN
207 208
                     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
209 210 211 212 213 214
                        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
215
                           DO l = ABS(l1 - l2), MIN(hybinp%lcutm1(itype), l1 + l2)
Matthias Redies's avatar
Matthias Redies committed
216 217 218
                              IF (MOD(l1 + l2 + l, 2) == 0) THEN
                                 hybdat%nindxp1(l, itype) = hybdat%nindxp1(l, itype) + 1
                                 n = hybdat%nindxp1(l, itype)
219 220 221 222 223 224
                                 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
225
                                                                          atoms, itype, hybdat%gridf)
Matthias Redies's avatar
Matthias Redies committed
226 227 228
                                 END DO
                              END IF
                           END DO
229 230
                        END DO
                     END DO
Matthias Redies's avatar
Matthias Redies committed
231 232
                  END IF
               END DO
233 234
            END DO
         END DO
235
         deallocate(basprod)
Matthias Redies's avatar
Matthias Redies committed
236
         CALL timestop("gen_bz and gen_wavf")
237

238
      ELSE IF (hybinp%l_hybrid) THEN ! hybdat%l_calhf is false
Matthias Redies's avatar
Matthias Redies committed
239 240 241

         !DO nk = n_start,kpts%nkpt,n_stride
         DO nk = 1, kpts%nkpt, 1
242 243
            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
244
         END DO
245

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

249
      ENDIF ! hybdat%l_calhf
250

Matthias Redies's avatar
Matthias Redies committed
251
   END SUBROUTINE hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
252

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