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

11
   SUBROUTINE hf_setup(mpbasis, hybrid, input, sym, kpts,  atoms, mpi, noco, cell, oneD, results, jsp, enpara, eig_id_hf, &
12
                       hybdat, l_real, vr0, eig_irr)
Matthias Redies's avatar
Matthias Redies committed
13 14 15
      USE m_types
      USE m_eig66_io
      USE m_util
Matthias Redies's avatar
Matthias Redies committed
16
      USE m_intgrf
Matthias Redies's avatar
Matthias Redies committed
17
      USE m_checkolap
18
      USE m_hybrid_core
Matthias Redies's avatar
Matthias Redies committed
19 20 21 22
      USE m_gen_wavf

      IMPLICIT NONE

Matthias Redies's avatar
Matthias Redies committed
23
      TYPE(t_mpbasis), INTENT(inout)   :: mpbasis
Matthias Redies's avatar
Matthias Redies committed
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
      TYPE(t_hybrid), INTENT(INOUT) :: hybrid
      TYPE(t_kpts), INTENT(IN)    :: kpts
      TYPE(t_atoms), INTENT(IN)    :: atoms
      TYPE(t_mpi), INTENT(IN)    :: mpi
      TYPE(t_noco), INTENT(IN)    :: noco
      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

      INTEGER, INTENT(IN)    :: jsp, eig_id_hf
      REAL, INTENT(IN)    :: vr0(:, :, :)
      LOGICAL, INTENT(IN)    :: l_real

      REAL, ALLOCATABLE, INTENT(OUT)   :: eig_irr(:, :)

      ! local type variables
      TYPE(t_lapw)             :: lapw
      TYPE(t_mat), ALLOCATABLE :: zmat(:)

      ! 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
49
      INTEGER :: nbasfcn, n_dim
50
      LOGICAL :: l_exist
Matthias Redies's avatar
Matthias Redies committed
51 52 53 54

      ! local arrays

      REAL, ALLOCATABLE :: basprod(:)
55
      INTEGER              :: degenerat(merge(input%neig*2,input%neig,noco%l_soc) + 1, kpts%nkpt)
Matthias Redies's avatar
Matthias Redies committed
56
      LOGICAL              :: skip_kpt(kpts%nkpt)
57

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

Matthias Redies's avatar
Matthias Redies committed
61
      skip_kpt = .FALSE.
Daniel Wortmann's avatar
Daniel Wortmann committed
62

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

67
         allocate(zmat(kpts%nkptf), stat=ok)
68
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation z_c')
69
         allocate(eig_irr(input%neig, kpts%nkpt), stat=ok)
70
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation eig_irr')
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 75
         eig_irr = 0
         hybdat%kveclo_eig = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
76

77 78
         INQUIRE(file ="z",exist= l_exist)
         IF(l_exist) THEN
79 80
            IF (l_real) OPEN(unit=993,file='z',form='unformatted',access='direct',recl=lapw%dim_nbasfcn()*input%neig*8)
            IF (.NOT.l_real) OPEN(unit=993,file='z',form='unformatted',access='direct',recl=lapw%dim_nbasfcn()*input%neig*16)
81 82
         END IF

Matthias Redies's avatar
Matthias Redies committed
83 84 85 86 87
         ! Reading the eig file
         DO nk = 1, kpts%nkpt
            nrec1 = kpts%nkpt*(jsp - 1) + nk
            CALL lapw%init(input, noco, kpts, atoms, sym, nk, cell, sym%zrfs)
            nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
88
            CALL zMat(nk)%init(l_real, nbasfcn, merge(input%neig*2,input%neig,noco%l_soc))
Matthias Redies's avatar
Matthias Redies committed
89
            CALL read_eig(eig_id_hf, nk, jsp, zmat=zMat(nk))
90

91 92 93
            IF(l_exist.AND.zmat(1)%l_real) THEN
               READ(993,rec=nk) zDebug_r(:,:)
               zMat(nk)%data_r = 0.0
94
               zMat(nk)%data_r(:nbasfcn,:input%neig) = zDebug_r(:nbasfcn,:input%neig)
95 96 97 98
            END IF
            IF(l_exist.AND..NOT.zmat(1)%l_real) THEN
               READ(993,rec=nk) zDebug_c(:,:)
               zMat(nk)%data_c = 0.0
99
               zMat(nk)%data_c(:nbasfcn,:input%neig) = zDebug_c(:nbasfcn,:input%neig)
100 101
            END IF

Matthias Redies's avatar
Matthias Redies committed
102 103 104
            eig_irr(:, nk) = results%eig(:, nk, jsp)
            hybrid%ne_eig(nk) = results%neig(nk, jsp)
         END DO
105 106 107

         IF(l_exist) CLOSE(993)

Matthias Redies's avatar
Matthias Redies committed
108 109 110
         !Allocate further space
         DO nk = kpts%nkpt + 1, kpts%nkptf
            nbasfcn = zMat(kpts%bkp(nk))%matsize1
111
            CALL zMat(nk)%init(l_real, nbasfcn, merge(input%neig*2,input%neig,noco%l_soc))
Matthias Redies's avatar
Matthias Redies committed
112 113 114 115 116 117 118 119 120 121 122 123 124
         END DO

         !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
125
         hybrid%nobd(:,jsp) = 0
Matthias Redies's avatar
Matthias Redies committed
126 127 128 129 130 131 132
         DO nk = 1, kpts%nkpt
            DO i = 1, hybrid%ne_eig(nk)
               DO j = i + 1, hybrid%ne_eig(nk)
                  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
133 134
            END DO

Matthias Redies's avatar
Matthias Redies committed
135 136 137
            DO i = 1, hybrid%ne_eig(nk)
               IF ((degenerat(i, nk) /= 1) .OR. (degenerat(i, nk) /= 0)) degenerat(i + 1:i + degenerat(i, nk) - 1, nk) = 0
            END DO
138

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

Matthias Redies's avatar
Matthias Redies committed
141 142 143 144 145 146 147 148
            hybrid%nbands(nk) = hybrid%bands1
            IF (hybrid%nbands(nk) > hybrid%ne_eig(nk)) THEN
               IF (mpi%irank == 0) THEN
                  WRITE (*, *) ' maximum for hybrid%nbands is', hybrid%ne_eig(nk)
                  WRITE (*, *) ' increase energy window to obtain enough eigenvalues'
                  WRITE (*, *) ' set hybrid%nbands equal to hybrid%ne_eig'
               END IF
               hybrid%nbands(nk) = hybrid%ne_eig(nk)
149 150
            END IF

Matthias Redies's avatar
Matthias Redies committed
151 152 153 154 155 156
            DO i = hybrid%nbands(nk) - 1, 1, -1
               IF ((degenerat(i, nk) >= 1) .AND. (degenerat(i, nk) + i - 1 /= hybrid%nbands(nk))) THEN
                  hybrid%nbands(nk) = i + degenerat(i, nk) - 1
                  EXIT
               END IF
            END DO
157

Matthias Redies's avatar
Matthias Redies committed
158
            DO i = 1, hybrid%ne_eig(nk)
159
               IF (results%w_iks(i, nk, jsp) > 0.0) hybrid%nobd(nk,jsp) = hybrid%nobd(nk,jsp) + 1
Matthias Redies's avatar
Matthias Redies committed
160
            END DO
161

162
            IF (hybrid%nobd(nk,jsp) > hybrid%nbands(nk)) THEN
Matthias Redies's avatar
Matthias Redies committed
163 164
               WRITE (*, *) 'k-point: ', nk
               WRITE (*, *) 'number of bands:          ', hybrid%nbands(nk)
165
               WRITE (*, *) 'number of occupied bands: ', hybrid%nobd(nk,jsp)
Matthias Redies's avatar
Matthias Redies committed
166
               CALL judft_warn("More occupied bands than total no of bands!?")
167
               hybrid%nbands(nk) = hybrid%nobd(nk,jsp)
Matthias Redies's avatar
Matthias Redies committed
168
            END IF
169
            PRINT *, "bands:", nk, hybrid%nobd(nk,jsp), hybrid%nbands(nk), hybrid%ne_eig(nk)
Matthias Redies's avatar
Matthias Redies committed
170
         END DO
171

Matthias Redies's avatar
Matthias Redies committed
172 173 174
         ! spread hybrid%nobd from IBZ to whole BZ
         DO nk = 1, kpts%nkptf
            i = kpts%bkp(nk)
175
            hybrid%nobd(nk,jsp) = hybrid%nobd(i,jsp)
Matthias Redies's avatar
Matthias Redies committed
176
         END DO
177

Matthias Redies's avatar
Matthias Redies committed
178
         ! generate eigenvectors z and MT coefficients from the previous iteration at all k-points
179
         CALL gen_wavf(kpts%nkpt, kpts, sym, atoms, enpara%el0(:, :, jsp), enpara%ello0(:, :, jsp), cell,  &
180
                       mpbasis, hybrid, vr0, hybdat, noco, oneD, mpi, input, jsp, zmat)
181

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

Matthias Redies's avatar
Matthias Redies committed
186
         ! check olap between core-basis/core-valence/basis-basis
187
         CALL checkolap(atoms, hybdat, mpbasis, hybrid, kpts%nkpt, kpts,  mpi, &
Matthias Redies's avatar
Matthias Redies committed
188 189 190 191 192
                        input, sym, noco, cell, lapw, jsp)

         ! set up pointer pntgpt

         ! setup dimension of pntgpt
193
         IF(ALLOCATED(hybdat%pntgptd)) DEALLOCATE(hybdat%pntgptd) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
194 195 196 197
         ALLOCATE (hybdat%pntgptd(3))
         hybdat%pntgptd = 0
         DO nk = 1, kpts%nkptf
            CALL lapw%init(input, noco, kpts, atoms, sym, nk, cell, sym%zrfs)
Matthias Redies's avatar
Matthias Redies committed
198
            do n_dim = 1,3
Matthias Redies's avatar
again  
Matthias Redies committed
199
               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
200
            end do
Matthias Redies's avatar
Matthias Redies committed
201 202
         END DO

203
         IF(ALLOCATED(hybdat%pntgpt)) DEALLOCATE(hybdat%pntgpt) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
204 205
         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)
206
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation pntgpt')
Matthias Redies's avatar
Matthias Redies committed
207 208 209 210
         hybdat%pntgpt = 0
         DO nk = 1, kpts%nkptf
            CALL lapw%init(input, noco, kpts, atoms, sym, nk, cell, sym%zrfs)
            DO i = 1, lapw%nv(jsp)
211
               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
212
            END DO
213
         END DO
Matthias Redies's avatar
Matthias Redies committed
214

215
         allocate(basprod(atoms%jmtd), stat=ok)
216
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation basprod')
Matthias Redies's avatar
Matthias Redies committed
217
         IF(ALLOCATED(hybdat%prodm)) DEALLOCATE(hybdat%prodm)
218
         allocate(hybdat%prodm(maxval(mpbasis%num_radbasfn), hybrid%max_indx_p_1, 0:maxval(hybrid%lcutm1), atoms%ntype), stat=ok)
219
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%prodm')
220

Matthias Redies's avatar
Matthias Redies committed
221
         call mpbasis%init(hybrid, atoms)
222

Matthias Redies's avatar
Matthias Redies committed
223 224
         basprod = 0; hybdat%prodm = 0; mpbasis%l1 = 0; mpbasis%l2 = 0
         mpbasis%n1 = 0; mpbasis%n2 = 0
225
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE(hybdat%nindxp1) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
226
         ALLOCATE (hybdat%nindxp1(0:maxval(hybrid%lcutm1), atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
227 228 229 230 231 232 233
         hybdat%nindxp1 = 0
         DO itype = 1, atoms%ntype
            ng = atoms%jri(itype)
            DO l2 = 0, MIN(atoms%lmax(itype), hybrid%lcutwf(itype))
               ll = l2
               DO l1 = 0, ll
                  IF (ABS(l1 - l2) <= hybrid%lcutm1(itype)) THEN
234 235
                     DO n2 = 1, mpbasis%num_radfun_per_l(l2, itype)
                        nn = mpbasis%num_radfun_per_l(l1, itype)
Matthias Redies's avatar
Matthias Redies committed
236 237 238 239 240 241 242 243 244 245
                        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)
                           DO l = ABS(l1 - l2), MIN(hybrid%lcutm1(itype), l1 + l2)
                              IF (MOD(l1 + l2 + l, 2) == 0) THEN
                                 hybdat%nindxp1(l, itype) = hybdat%nindxp1(l, itype) + 1
                                 n = hybdat%nindxp1(l, itype)
Matthias Redies's avatar
Matthias Redies committed
246 247 248 249
                                 mpbasis%l1(n,l,itype) = l1
                                 mpbasis%l2(n,l,itype) = l2
                                 mpbasis%n1(n,l,itype) = n1
                                 mpbasis%n2(n,l,itype) = n2
250
                                 DO i = 1, mpbasis%num_radbasfn(l, itype)
Matthias Redies's avatar
Matthias Redies committed
251 252
                                    hybdat%prodm(i, n, l, itype) = intgrf(basprod(:ng)*mpbasis%radbasfn_mt(:ng, i, l, itype), &
                                                                          atoms, itype, hybdat%gridf)
Matthias Redies's avatar
Matthias Redies committed
253 254 255
                                 END DO
                              END IF
                           END DO
256 257
                        END DO
                     END DO
Matthias Redies's avatar
Matthias Redies committed
258 259
                  END IF
               END DO
260 261
            END DO
         END DO
262
         deallocate(basprod)
Matthias Redies's avatar
Matthias Redies committed
263
         CALL timestop("gen_bz and gen_wavf")
264

Matthias Redies's avatar
Matthias Redies committed
265 266 267 268 269
      ELSE IF (hybrid%l_hybrid) THEN ! hybrid%l_calhf is false

         !DO nk = n_start,kpts%nkpt,n_stride
         DO nk = 1, kpts%nkpt, 1
            hybrid%ne_eig(nk) = results%neig(nk, jsp)
270
            hybrid%nobd(nk,jsp) = COUNT(results%w_iks(:hybrid%ne_eig(nk), nk, jsp) > 0.0)
Matthias Redies's avatar
Matthias Redies committed
271
         END DO
272

273
         hybrid%maxlmindx = MAXVAL([(SUM([(mpbasis%num_radfun_per_l(l, itype)*(2*l + 1), l=0, atoms%lmax(itype))]), itype=1, atoms%ntype)])
274
         hybrid%nbands = MIN(hybrid%bands1, input%neig)
275

Matthias Redies's avatar
Matthias Redies committed
276
      ENDIF ! hybrid%l_calhf
277

Matthias Redies's avatar
Matthias Redies committed
278
   END SUBROUTINE hf_setup
Daniel Wortmann's avatar
Daniel Wortmann committed
279

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