hf_setup.F90 11.4 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
      USE m_types
15
      USE m_constants
Matthias Redies's avatar
Matthias Redies committed
16 17
      USE m_eig66_io
      USE m_util
Matthias Redies's avatar
Matthias Redies committed
18
      USE m_intgrf
Matthias Redies's avatar
Matthias Redies committed
19
      USE m_checkolap
Matthias Redies's avatar
Matthias Redies committed
20
      USE m_hybrid_core
Matthias Redies's avatar
Matthias Redies committed
21
      USE m_gen_wavf
Daniel Wortmann's avatar
Daniel Wortmann committed
22
      use m_types_hybdat
Matthias Redies's avatar
Matthias Redies committed
23 24 25

      IMPLICIT NONE

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

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

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

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

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

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

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

78

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

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

         !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
Matthias Redies's avatar
Matthias Redies committed
97 98

         call timestart("degenerate treatment")
Matthias Redies's avatar
Matthias Redies committed
99
         IF (mpi%irank == 0) THEN
100 101
            WRITE (oUnit, *)
            WRITE (oUnit, '(A)') "   k-point      |   number of occupied bands  |   maximal number of bands"
Matthias Redies's avatar
Matthias Redies committed
102 103
         END IF
         degenerat = 1
Matthias Redies's avatar
Matthias Redies committed
104
         hybdat%nobd(:,jsp) = 0         
Matthias Redies's avatar
Matthias Redies committed
105
         DO nk = 1, kpts%nkpt
106 107
            DO i = 1, hybdat%ne_eig(nk)
               DO j = i + 1, hybdat%ne_eig(nk)
Matthias Redies's avatar
Matthias Redies committed
108 109 110 111
                  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
112 113
            END DO

114
            DO i = 1, hybdat%ne_eig(nk)
Matthias Redies's avatar
Matthias Redies committed
115 116
               IF ((degenerat(i, nk) /= 1) .OR. (degenerat(i, nk) /= 0)) degenerat(i + 1:i + degenerat(i, nk) - 1, nk) = 0
            END DO
117

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

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

130 131 132
            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
133 134 135
                  EXIT
               END IF
            END DO
136

137 138
            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
139
            END DO
140

141
            IF (hybdat%nobd(nk,jsp) > hybdat%nbands(nk)) THEN
Matthias Redies's avatar
Matthias Redies committed
142
               WRITE (*, *) 'k-point: ', nk
143 144
               WRITE (*, *) 'number of bands:          ', hybdat%nbands(nk)
               WRITE (*, *) 'number of occupied bands: ', hybdat%nobd(nk,jsp)
Matthias Redies's avatar
Matthias Redies committed
145
               CALL judft_warn("More occupied bands than total no of bands!?")
146
               hybdat%nbands(nk) = hybdat%nobd(nk,jsp)
Matthias Redies's avatar
Matthias Redies committed
147 148
            END IF
         END DO
Matthias Redies's avatar
Matthias Redies committed
149
         call timestop("degenerate treatment")
150

151
         ! spread hybdat%nobd from IBZ to whole BZ
Matthias Redies's avatar
Matthias Redies committed
152 153
         DO nk = 1, kpts%nkptf
            i = kpts%bkp(nk)
Matthias Redies's avatar
Matthias Redies committed
154
            hybdat%nbands(nk) = hybdat%nbands(i)
155
            hybdat%nobd(nk,jsp) = hybdat%nobd(i,jsp)
Matthias Redies's avatar
Matthias Redies committed
156
         END DO
157

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

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

Matthias Redies's avatar
Matthias Redies committed
166
         ! check olap between core-basis/core-valence/basis-basis
Matthias Redies's avatar
Matthias Redies committed
167
         CALL checkolap(atoms, hybdat, mpdata, hybinp, kpts%nkpt, kpts,  mpi, &
Matthias Redies's avatar
Matthias Redies committed
168
                        input, sym, noco, nococonv,oneD,cell, lapw, jsp)
Matthias Redies's avatar
Matthias Redies committed
169 170 171 172

         ! set up pointer pntgpt

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

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

Matthias Redies's avatar
Matthias Redies committed
195
         allocate(basprod(atoms%jmtd), stat=ok, source=0.0)
196
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation basprod')
Matthias Redies's avatar
Matthias Redies committed
197
         IF(ALLOCATED(hybdat%prodm)) DEALLOCATE(hybdat%prodm)
198
         allocate(hybdat%prodm(maxval(mpdata%num_radbasfn), mpdata%max_indx_p_1, 0:maxval(hybinp%lcutm1), atoms%ntype), stat=ok)
199
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%prodm')
200

Matthias Redies's avatar
Matthias Redies committed
201
         hybdat%prodm = 0; mpdata%l1 = 0; mpdata%l2 = 0
202
         mpdata%n1 = 0; mpdata%n2 = 0
203
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE(hybdat%nindxp1) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
204
         ALLOCATE (hybdat%nindxp1(0:maxval(hybinp%lcutm1), atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
205
         hybdat%nindxp1 = 0
Matthias Redies's avatar
Matthias Redies committed
206 207 208 209

         !$OMP PARALLEL DO default(none) schedule(dynamic)&
         !$OMP private(itype, ng, l2, l1, n1, l, nn, n, basprod) &
         !$OMP shared(atoms, hybinp, mpdata, hybdat)
Matthias Redies's avatar
Matthias Redies committed
210 211
         DO itype = 1, atoms%ntype
            ng = atoms%jri(itype)
Matthias Redies's avatar
Matthias Redies committed
212
            DO l2 = 0, MIN(atoms%lmax(itype), hybinp%lcutwf(itype))
Matthias Redies's avatar
Matthias Redies committed
213
               DO l1 = 0, l2
Matthias Redies's avatar
Matthias Redies committed
214
                  IF (ABS(l1 - l2) <= hybinp%lcutm1(itype)) THEN
215 216
                     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
217 218 219 220 221 222
                        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
223
                           DO l = ABS(l1 - l2), MIN(hybinp%lcutm1(itype), l1 + l2)
Matthias Redies's avatar
Matthias Redies committed
224 225 226
                              IF (MOD(l1 + l2 + l, 2) == 0) THEN
                                 hybdat%nindxp1(l, itype) = hybdat%nindxp1(l, itype) + 1
                                 n = hybdat%nindxp1(l, itype)
227 228 229 230 231 232
                                 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
233
                                                                          atoms, itype, hybdat%gridf)
Matthias Redies's avatar
Matthias Redies committed
234 235 236
                                 END DO
                              END IF
                           END DO
237 238
                        END DO
                     END DO
Matthias Redies's avatar
Matthias Redies committed
239 240
                  END IF
               END DO
241 242
            END DO
         END DO
Matthias Redies's avatar
Matthias Redies committed
243
         !$OMP END PARALLEL DO
244
         deallocate(basprod)
Matthias Redies's avatar
Matthias Redies committed
245
         CALL timestop("gen_bz and gen_wavf")
246

247
      ELSE IF (hybinp%l_hybrid) THEN ! hybdat%l_calhf is false
Matthias Redies's avatar
Matthias Redies committed
248 249 250

         !DO nk = n_start,kpts%nkpt,n_stride
         DO nk = 1, kpts%nkpt, 1
251 252
            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
253
         END DO
254

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

258
      ENDIF ! hybdat%l_calhf
259

Matthias Redies's avatar
Matthias Redies committed
260
   END SUBROUTINE hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
261

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