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

Matthias Redies's avatar
Matthias Redies committed
7

Daniel Wortmann's avatar
Daniel Wortmann committed
8
MODULE m_hf_setup
9

Daniel Wortmann's avatar
Daniel Wortmann committed
10
CONTAINS
11

Matthias Redies's avatar
Matthias Redies committed
12
   SUBROUTINE hf_setup(mpdata, fi, mpi,nococonv, 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

Matthias Redies's avatar
Matthias Redies committed
26
      type(t_fleurinput), intent(in)    :: fi
27
      TYPE(t_mpdata), INTENT(inout)   :: mpdata
Matthias Redies's avatar
Matthias Redies committed
28
      TYPE(t_mpi), INTENT(IN)    :: mpi
29
      TYPE(t_nococonv), INTENT(IN)    :: nococonv
Matthias Redies's avatar
Matthias Redies committed
30 31 32 33
      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
34
      INTEGER, INTENT(IN)    :: jsp
Matthias Redies's avatar
Matthias Redies committed
35 36 37
      REAL, INTENT(IN)    :: vr0(:, :, :)
      LOGICAL, INTENT(IN)    :: l_real

Matthias Redies's avatar
Matthias Redies committed
38
      REAL, ALLOCATABLE, INTENT(INOUT)   :: eig_irr(:, :)
Matthias Redies's avatar
Matthias Redies committed
39 40 41 42 43 44

      ! 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
45
      INTEGER :: nbasfcn, n_dim
46
      LOGICAL :: l_exist
Matthias Redies's avatar
Matthias Redies committed
47 48 49 50

      ! local arrays

      REAL, ALLOCATABLE :: basprod(:)
Matthias Redies's avatar
Matthias Redies committed
51
      INTEGER              :: degenerat(merge(input%neig*2,input%neig,noco%l_soc) + 1, fi%kpts%nkpt)
52

53 54
      REAL :: zDebug_r(lapw_dim_nbasfcn,input%neig)
      COMPLEX :: zDebug_c(lapw_dim_nbasfcn,input%neig)
55

56
      IF (hybdat%l_calhf) THEN
Matthias Redies's avatar
Matthias Redies committed
57
         ! Preparations for HF and hybinp functional calculation
Matthias Redies's avatar
Matthias Redies committed
58
         CALL timestart("gen_bz and gen_wavf")
59

Matthias Redies's avatar
Matthias Redies committed
60
         if(.not. allocated(eig_irr)) then 
Matthias Redies's avatar
Matthias Redies committed
61
            allocate(eig_irr(input%neig, fi%kpts%nkpt), stat=ok)
Matthias Redies's avatar
Matthias Redies committed
62 63 64
            IF (ok /= 0) call judft_error('eigen_hf: failure allocation eig_irr')
         endif
         eig_irr = 0.0
Matthias Redies's avatar
Matthias Redies committed
65

Matthias Redies's avatar
Matthias Redies committed
66
         if(allocated(hybdat%kveclo_eig)) deallocate(hybdat%kveclo_eig)
Matthias Redies's avatar
Matthias Redies committed
67
         allocate(hybdat%kveclo_eig(fi%atoms%nlotot, fi%kpts%nkpt), stat=ok)
68
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%kveclo_eig')
Matthias Redies's avatar
Matthias Redies committed
69
         hybdat%kveclo_eig = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
70

71

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

Matthias Redies's avatar
Matthias Redies committed
79
            eig_irr(:, nk) = results%eig(:, nk, jsp)
80
            hybdat%ne_eig(nk) = results%neig(nk, jsp)
Matthias Redies's avatar
Matthias Redies committed
81
         END DO
Matthias Redies's avatar
Matthias Redies committed
82
         call timestop("eig stuff")
Matthias Redies's avatar
Matthias Redies committed
83 84 85 86 87 88 89

         !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
90 91

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

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

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

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

123 124 125
            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
126 127 128
                  EXIT
               END IF
            END DO
129

130 131
            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
132
            END DO
133

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

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

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

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

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

         ! set up pointer pntgpt

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

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

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

Matthias Redies's avatar
Matthias Redies committed
194
         hybdat%prodm = 0; mpdata%l1 = 0; mpdata%l2 = 0
195
         mpdata%n1 = 0; mpdata%n2 = 0
196
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE(hybdat%nindxp1) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
197
         ALLOCATE (hybdat%nindxp1(0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
198
         hybdat%nindxp1 = 0
Matthias Redies's avatar
Matthias Redies committed
199 200 201

         !$OMP PARALLEL DO default(none) schedule(dynamic)&
         !$OMP private(itype, ng, l2, l1, n1, l, nn, n, basprod) &
Matthias Redies's avatar
Matthias Redies committed
202 203 204 205
         !$OMP shared(fi, mpdata, hybdat)
         DO itype = 1, fi%atoms%ntype
            ng = fi%atoms%jri(itype)
            DO l2 = 0, MIN(fi%atoms%lmax(itype), fi%hybinp%lcutwf(itype))
Matthias Redies's avatar
Matthias Redies committed
206
               DO l1 = 0, l2
Matthias Redies's avatar
Matthias Redies committed
207
                  IF (ABS(l1 - l2) <= fi%hybinp%lcutm1(itype)) THEN
208 209
                     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
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) + &
Matthias Redies's avatar
Matthias Redies committed
215 216
                                           hybdat%bas2(:ng, n1, l1, itype)*hybdat%bas2(:ng, n2, l2, itype))/fi%atoms%rmsh(:ng, itype)
                           DO l = ABS(l1 - l2), MIN(fi%hybinp%lcutm1(itype), l1 + l2)
Matthias Redies's avatar
Matthias Redies committed
217 218 219
                              IF (MOD(l1 + l2 + l, 2) == 0) THEN
                                 hybdat%nindxp1(l, itype) = hybdat%nindxp1(l, itype) + 1
                                 n = hybdat%nindxp1(l, itype)
220 221 222 223 224 225
                                 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
226
                                                                          fi%atoms, itype, hybdat%gridf)
Matthias Redies's avatar
Matthias Redies committed
227 228 229
                                 END DO
                              END IF
                           END DO
230 231
                        END DO
                     END DO
Matthias Redies's avatar
Matthias Redies committed
232 233
                  END IF
               END DO
234 235
            END DO
         END DO
Matthias Redies's avatar
Matthias Redies committed
236
         !$OMP END PARALLEL DO
237
         deallocate(basprod)
Matthias Redies's avatar
Matthias Redies committed
238
         CALL timestop("gen_bz and gen_wavf")
239

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

Matthias Redies's avatar
Matthias Redies committed
242 243
         !DO nk = n_start,fi%kpts%nkpt,n_stride
         DO nk = 1, fi%kpts%nkpt, 1
244 245
            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
246
         END DO
247

Matthias Redies's avatar
Matthias Redies committed
248 249
         hybdat%maxlmindx = MAXVAL([(SUM([(mpdata%num_radfun_per_l(l, itype)*(2*l + 1), l=0, fi%atoms%lmax(itype))]), itype=1, fi%atoms%ntype)])
         hybdat%nbands = MIN(fi%hybinp%bands1, input%neig)
250

251
      ENDIF ! hybdat%l_calhf
252

Matthias Redies's avatar
Matthias Redies committed
253
   END SUBROUTINE hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
254

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