hf_setup.F90 12.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, cell,&
                       oneD, results, jsp, enpara, 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
Matthias Redies's avatar
Matthias Redies committed
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
25
      TYPE(t_hybinp), INTENT(IN) :: hybinp
Matthias Redies's avatar
Matthias Redies committed
26 27 28 29 30 31 32 33 34 35 36 37
      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

Matthias Redies's avatar
Matthias Redies committed
38
      INTEGER, INTENT(IN)    :: jsp
Matthias Redies's avatar
Matthias Redies committed
39 40 41
      REAL, INTENT(IN)    :: vr0(:, :, :)
      LOGICAL, INTENT(IN)    :: l_real

Matthias Redies's avatar
Matthias Redies committed
42
      REAL, ALLOCATABLE, INTENT(INOUT)   :: eig_irr(:, :)
Matthias Redies's avatar
Matthias Redies committed
43 44 45 46 47 48 49

      ! 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

64
      IF (hybdat%l_calhf) THEN
Matthias Redies's avatar
Matthias Redies committed
65
         ! Preparations for HF and hybinp functional calculation
Matthias Redies's avatar
Matthias Redies committed
66
         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 71

         if(.not. allocated(eig_irr)) allocate(eig_irr(input%neig, kpts%nkpt), stat=ok)
72
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation eig_irr')
73
         eig_irr = 0
Matthias Redies's avatar
Matthias Redies committed
74

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

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

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

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

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

         IF(l_exist) CLOSE(993)

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

138
            DO i = 1, hybdat%ne_eig(nk)
Matthias Redies's avatar
Matthias Redies committed
139 140
               IF ((degenerat(i, nk) /= 1) .OR. (degenerat(i, nk) /= 0)) degenerat(i + 1:i + degenerat(i, nk) - 1, nk) = 0
            END DO
141

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

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

154 155 156
            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
157 158 159
                  EXIT
               END IF
            END DO
160

161 162
            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
163
            END DO
164

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

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

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

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

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

         ! set up pointer pntgpt

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

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

218
         allocate(basprod(atoms%jmtd), stat=ok)
219
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation basprod')
Matthias Redies's avatar
Matthias Redies committed
220
         IF(ALLOCATED(hybdat%prodm)) DEALLOCATE(hybdat%prodm)
221
         allocate(hybdat%prodm(maxval(mpdata%num_radbasfn), mpdata%max_indx_p_1, 0:maxval(hybinp%lcutm1), atoms%ntype), stat=ok)
222
         IF (ok /= 0) call judft_error('eigen_hf: failure allocation hybdat%prodm')
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(hybinp%lcutm1), atoms%ntype))
Matthias Redies's avatar
Matthias Redies committed
228 229 230
         hybdat%nindxp1 = 0
         DO itype = 1, atoms%ntype
            ng = atoms%jri(itype)
Matthias Redies's avatar
Matthias Redies committed
231
            DO l2 = 0, MIN(atoms%lmax(itype), hybinp%lcutwf(itype))
Matthias Redies's avatar
Matthias Redies committed
232 233
               ll = l2
               DO l1 = 0, ll
Matthias Redies's avatar
Matthias Redies committed
234
                  IF (ABS(l1 - l2) <= hybinp%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
                        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
243
                           DO l = ABS(l1 - l2), MIN(hybinp%lcutm1(itype), l1 + l2)
Matthias Redies's avatar
Matthias Redies committed
244 245 246
                              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

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

         !DO nk = n_start,kpts%nkpt,n_stride
         DO nk = 1, kpts%nkpt, 1
270 271
            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
272
         END DO
273

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

277
      ENDIF ! hybdat%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