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(mpdata, 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
      USE m_gen_wavf
Daniel Wortmann's avatar
Daniel Wortmann committed
20
      use m_types_hybdat
Matthias Redies's avatar
Matthias Redies committed
21 22 23

      IMPLICIT NONE

24
      TYPE(t_mpdata), INTENT(inout)   :: mpdata
Matthias Redies's avatar
Matthias Redies committed
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
      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
50
      INTEGER :: nbasfcn, n_dim
51
      LOGICAL :: l_exist
Matthias Redies's avatar
Matthias Redies committed
52 53 54 55

      ! local arrays

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

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

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

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

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

78 79
         INQUIRE(file ="z",exist= l_exist)
         IF(l_exist) THEN
80 81
            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)
82 83
         END IF

Matthias Redies's avatar
Matthias Redies committed
84 85 86 87 88
         ! 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)
89
            CALL zMat(nk)%init(l_real, nbasfcn, merge(input%neig*2,input%neig,noco%l_soc))
Matthias Redies's avatar
Matthias Redies committed
90
            CALL read_eig(eig_id_hf, nk, jsp, zmat=zMat(nk))
91

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

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

         IF(l_exist) CLOSE(993)

Matthias Redies's avatar
Matthias Redies committed
109 110 111
         !Allocate further space
         DO nk = kpts%nkpt + 1, kpts%nkptf
            nbasfcn = zMat(kpts%bkp(nk))%matsize1
112
            CALL zMat(nk)%init(l_real, nbasfcn, merge(input%neig*2,input%neig,noco%l_soc))
Matthias Redies's avatar
Matthias Redies committed
113 114 115 116 117 118 119 120 121 122 123 124 125
         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
126
         hybrid%nobd(:,jsp) = 0
Matthias Redies's avatar
Matthias Redies committed
127 128 129 130 131 132 133
         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
134 135
            END DO

Matthias Redies's avatar
Matthias Redies committed
136 137 138
            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
139

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

Matthias Redies's avatar
Matthias Redies committed
142 143 144 145 146 147 148 149
            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)
150 151
            END IF

Matthias Redies's avatar
Matthias Redies committed
152 153 154 155 156 157
            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
158

Matthias Redies's avatar
Matthias Redies committed
159
            DO i = 1, hybrid%ne_eig(nk)
160
               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
161
            END DO
162

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

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

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

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

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

         ! set up pointer pntgpt

         ! setup dimension of pntgpt
194
         IF(ALLOCATED(hybdat%pntgptd)) DEALLOCATE(hybdat%pntgptd) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
195 196 197 198
         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
199
            do n_dim = 1,3
Matthias Redies's avatar
again  
Matthias Redies committed
200
               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
201
            end do
Matthias Redies's avatar
Matthias Redies committed
202 203
         END DO

204
         IF(ALLOCATED(hybdat%pntgpt)) DEALLOCATE(hybdat%pntgpt) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
205 206
         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)
207
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation pntgpt')
Matthias Redies's avatar
Matthias Redies committed
208 209 210 211
         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)
212
               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
213
            END DO
214
         END DO
Matthias Redies's avatar
Matthias Redies committed
215

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

222
         call mpdata%init(hybrid, atoms)
223

224 225
         basprod = 0; hybdat%prodm = 0; mpdata%l1 = 0; mpdata%l2 = 0
         mpdata%n1 = 0; mpdata%n2 = 0
226
         IF(ALLOCATED(hybdat%nindxp1)) DEALLOCATE(hybdat%nindxp1) ! for spinpolarized systems
Matthias Redies's avatar
Matthias Redies committed
227
         ALLOCATE (hybdat%nindxp1(0:maxval(hybrid%lcutm1), atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
228 229 230 231 232 233 234
         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
235 236
                     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
237 238 239 240 241 242 243 244 245 246
                        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)
247 248 249 250 251 252
                                 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
253
                                                                          atoms, itype, hybdat%gridf)
Matthias Redies's avatar
Matthias Redies committed
254 255 256
                                 END DO
                              END IF
                           END DO
257 258
                        END DO
                     END DO
Matthias Redies's avatar
Matthias Redies committed
259 260
                  END IF
               END DO
261 262
            END DO
         END DO
263
         deallocate(basprod)
Matthias Redies's avatar
Matthias Redies committed
264
         CALL timestop("gen_bz and gen_wavf")
265

Matthias Redies's avatar
Matthias Redies committed
266 267 268 269 270
      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)
271
            hybrid%nobd(nk,jsp) = COUNT(results%w_iks(:hybrid%ne_eig(nk), nk, jsp) > 0.0)
Matthias Redies's avatar
Matthias Redies committed
272
         END DO
273

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

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

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

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