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

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

Matthias Redies's avatar
Matthias Redies committed
66 67 68 69 70
         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
71

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

77

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

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

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

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

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

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

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

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

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

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

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

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

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

         ! set up pointer pntgpt

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

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

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

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

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

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

         !DO nk = n_start,kpts%nkpt,n_stride
         DO nk = 1, kpts%nkpt, 1
247 248
            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
249
         END DO
250

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

254
      ENDIF ! hybdat%l_calhf
255

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

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