hf_setup.F90 11.3 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(fi%input%neig*2,fi%input%neig,fi%noco%l_soc) + 1, fi%kpts%nkpt)
52

Matthias Redies's avatar
Matthias Redies committed
53 54 55 56
      REAL :: zDebug_r(lapw_dim_nbasfcn,fi%input%neig)
      COMPLEX :: zDebug_c(lapw_dim_nbasfcn,fi%input%neig)

      call hybdat%set_states(fi, results, jsp)
57

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

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

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

73

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

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

         !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
92 93

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

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

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

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

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

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

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

146
         ! spread hybdat%nobd from IBZ to whole BZ
Matthias Redies's avatar
Matthias Redies committed
147 148
         DO nk = 1, fi%kpts%nkptf
            i = fi%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
Matthias Redies's avatar
Matthias Redies committed
154 155
         CALL gen_wavf(fi%kpts%nkpt, fi%kpts, fi%sym, fi%atoms, enpara%el0(:, :, jsp), enpara%ello0(:, :, jsp), fi%cell,  &
                       mpdata, fi%hybinp, vr0, hybdat, fi%noco, nococonv,fi%oneD, mpi, fi%input, jsp)
156

Matthias Redies's avatar
Matthias Redies committed
157
         ! generate core wave functions (-> core1/2(jmtd,hybdat%nindxc,0:lmaxc,ntype) )
Matthias Redies's avatar
Matthias Redies committed
158
         CALL corewf(fi%atoms, jsp, fi%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(fi%atoms, hybdat, mpdata, fi%hybinp, fi%kpts%nkpt, fi%kpts,  mpi, &
Matthias Redies's avatar
Matthias Redies committed
163
                        fi%input, fi%sym, fi%noco, nococonv,fi%oneD,fi%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
         ALLOCATE (hybdat%pntgptd(3))
         hybdat%pntgptd = 0
Matthias Redies's avatar
Matthias Redies committed
171
         DO nk = 1, fi%kpts%nkptf
Matthias Redies's avatar
Matthias Redies committed
172
            CALL lapw%init(fi%input, fi%noco, nococonv,fi%kpts, fi%atoms, fi%sym, nk, fi%cell, fi%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
         ALLOCATE (hybdat%pntgpt(-hybdat%pntgptd(1):hybdat%pntgptd(1), -hybdat%pntgptd(2):hybdat%pntgptd(2), &
Matthias Redies's avatar
Matthias Redies committed
180
                                 -hybdat%pntgptd(3):hybdat%pntgptd(3), fi%kpts%nkptf), stat=ok)
181
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation pntgpt')
Matthias Redies's avatar
Matthias Redies committed
182
         hybdat%pntgpt = 0
Matthias Redies's avatar
Matthias Redies committed
183
         DO nk = 1, fi%kpts%nkptf
Matthias Redies's avatar
Matthias Redies committed
184
            CALL lapw%init(fi%input, fi%noco, nococonv,fi%kpts, fi%atoms, fi%sym, nk, fi%cell, fi%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

Matthias Redies's avatar
Matthias Redies committed
190
         allocate(basprod(fi%atoms%jmtd), stat=ok, source=0.0)
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)
Matthias Redies's avatar
Matthias Redies committed
193
         allocate(hybdat%prodm(maxval(mpdata%num_radbasfn), mpdata%max_indx_p_1, 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype), stat=ok)
194
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%prodm')
195

Matthias Redies's avatar
Matthias Redies committed
196
         hybdat%prodm = 0; mpdata%l1 = 0; mpdata%l2 = 0
197
         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(fi%hybinp%lcutm1), fi%atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
200
         hybdat%nindxp1 = 0
Matthias Redies's avatar
Matthias Redies committed
201 202 203

         !$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
204 205 206 207
         !$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
208
               DO l1 = 0, l2
Matthias Redies's avatar
Matthias Redies committed
209
                  IF (ABS(l1 - l2) <= fi%hybinp%lcutm1(itype)) THEN
210 211
                     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
212 213 214 215 216
                        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
217 218
                                           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
219 220 221
                              IF (MOD(l1 + l2 + l, 2) == 0) THEN
                                 hybdat%nindxp1(l, itype) = hybdat%nindxp1(l, itype) + 1
                                 n = hybdat%nindxp1(l, itype)
222 223 224 225 226 227
                                 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
228
                                                                          fi%atoms, itype, hybdat%gridf)
Matthias Redies's avatar
Matthias Redies committed
229 230 231
                                 END DO
                              END IF
                           END DO
232 233
                        END DO
                     END DO
Matthias Redies's avatar
Matthias Redies committed
234 235
                  END IF
               END DO
236 237
            END DO
         END DO
Matthias Redies's avatar
Matthias Redies committed
238
         !$OMP END PARALLEL DO
239
         deallocate(basprod)
Matthias Redies's avatar
Matthias Redies committed
240
         CALL timestop("gen_bz and gen_wavf")
241

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

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

Matthias Redies's avatar
Matthias Redies committed
250
         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)])
Matthias Redies's avatar
Matthias Redies committed
251
         hybdat%nbands = MIN(fi%hybinp%bands1, fi%input%neig)
252

253
      ENDIF ! hybdat%l_calhf
254

Matthias Redies's avatar
Matthias Redies committed
255
   END SUBROUTINE hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
256

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