metagga.F90 14.8 KB
Newer Older
Matthias Redies's avatar
Matthias Redies committed
1 2 3 4 5
!--------------------------------------------------------------------------------
! Copyright (c) 2018 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
MODULE m_metagga
Matthias Redies's avatar
Matthias Redies committed
7
   PUBLIC  :: calc_EnergyDen
Matthias Redies's avatar
Matthias Redies committed
8 9 10
   PRIVATE :: calc_EnergyDen_auxillary_weights, &
              calc_kinEnergyDen_pw, &
              calc_kinEnergyDen_mt
Matthias Redies's avatar
Matthias Redies committed
11

Matthias Redies's avatar
Matthias Redies committed
12 13 14
   type t_RS_potden
      REAL, ALLOCATABLE :: is(:,:), mt(:,:)
   end type t_RS_potden
15

16 17 18 19 20 21 22 23
   TYPE t_kinED
      logical             :: set=.false.
      real, allocatable   :: is(:,:)   ! (nsp*jmtd, jspins)
      real, allocatable   :: mt(:,:,:) ! (nsp*jmtd, jspins, local num of types)
   contains
      procedure           :: alloc_mt => kED_alloc_mt
   END TYPE t_kinED

24
CONTAINS
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

  SUBROUTINE kED_alloc_mt(kED,nsp_x_jmtd, jspins, n_start, n_types, n_stride)
    IMPLICIT NONE
    class(t_kinED), intent(inout)   :: kED
    integer, intent(in)            :: nsp_x_jmtd, jspins, n_start, n_types, n_stride
    integer                        :: cnt, n

    if(.not. allocated(kED%mt)) then
      cnt = 0
      do n = n_start,n_types,n_stride
        cnt = cnt + 1
      enddo
      allocate(kED%mt(nsp_x_jmtd, jspins, cnt), source=0.0)
    endif
  END SUBROUTINE kED_alloc_mt


42
   SUBROUTINE calc_kinEnergyDen_pw(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS)
43
      USE m_juDFT_stop
Matthias Redies's avatar
Matthias Redies committed
44
      !use m_cdngen
Matthias Redies's avatar
Matthias Redies committed
45 46 47
      IMPLICIT NONE
      REAL, INTENT(in)                 :: den_RS(:,:), EnergyDen_RS(:,:), vTot_RS(:,:)
      REAL, INTENT(inout), allocatable :: kinEnergyDen_RS(:,:)
48
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
49
      REAL, PARAMETER                  :: eps = 1e-15
Matthias Redies's avatar
Matthias Redies committed
50 51

      kinEnergyDen_RS = EnergyDen_RS - vTot_RS * den_RS
52 53 54 55 56 57
#else
      CALL juDFT_error("MetaGGA require LibXC",hint="compile Fleur with LibXC (e.g. by giving '-external libxc' to ./configure")
#endif
   END SUBROUTINE calc_kinEnergyDen_pw

   SUBROUTINE calc_kinEnergyDen_mt(EnergyDen_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, &
Matthias Redies's avatar
Matthias Redies committed
58
                                   kinEnergyDen_RS)
59 60 61 62
      USE m_juDFT_stop
      USE m_juDFT_string
      implicit none
      REAL, INTENT(in)                 :: EnergyDen_RS(:,:), vTot_rs(:,:), vTot0_rs(:,:), core_den_rs(:,:), val_den_rs(:,:)
63
      REAL, INTENT(inout)              :: kinEnergyDen_RS(:,:)
64 65 66

#ifdef CPP_LIBXC
      kinEnergyDen_RS = EnergyDen_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
Matthias Redies's avatar
Matthias Redies committed
67
#else
Matthias Redies's avatar
Matthias Redies committed
68
      CALL juDFT_error("MetaGGA require LibXC",hint="compile Fleur with LibXC (e.g. by giving '-external libxc' to ./configure")
Matthias Redies's avatar
Matthias Redies committed
69
#endif
70
   END SUBROUTINE calc_kinEnergyDen_mt
Matthias Redies's avatar
Matthias Redies committed
71

Matthias Redies's avatar
Matthias Redies committed
72

Matthias Redies's avatar
Matthias Redies committed
73
   SUBROUTINE calc_EnergyDen(eig_id, mpi, kpts, noco, input, banddos, cell, atoms, enpara, stars, &
74
         vacuum,  sphhar, sym, gfinp, hub1inp, vTot, oneD, results, EnergyDen)
Matthias Redies's avatar
Matthias Redies committed
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
      ! calculates the energy density
      ! EnergyDen = \sum_i n_i(r) \varepsilon_i
      ! where n_i(r) is the one-particle density
      ! and \varepsilon_i are the eigenenergies

      USE m_types_setup
      USE m_types_potden
      USE m_types_kpts
      USE m_types_mpi
      USE m_types_enpara
      USE m_types_misc
      USE m_types_regionCharges
      USE m_types_dos
      USE m_types_cdnval
      USE m_cdnval

      IMPLICIT NONE

      INTEGER,           INTENT(in)           :: eig_id
      TYPE(t_mpi),       INTENT(in)           :: mpi
      TYPE(t_kpts),      INTENT(in)           :: kpts
      TYPE(t_noco),      INTENT(in)           :: noco
      TYPE(t_input),     INTENT(in)           :: input
      TYPE(t_banddos),   INTENT(in)           :: banddos
      TYPE(t_cell),      INTENT(in)           :: cell
      TYPE(t_atoms),     INTENT(in)           :: atoms
      TYPE(t_enpara),    INTENT(in)           :: enpara
      TYPE(t_stars),     INTENT(in)           :: stars
      TYPE(t_vacuum),    INTENT(in)           :: vacuum
104

Matthias Redies's avatar
Matthias Redies committed
105 106
      TYPE(t_sphhar),    INTENT(in)           :: sphhar
      TYPE(t_sym),       INTENT(in)           :: sym
107 108
      TYPE(t_gfinp),     INTENT(in)           :: gfinp
      TYPE(t_hub1inp),   INTENT(in)           :: hub1inp
Matthias Redies's avatar
Matthias Redies committed
109 110 111 112 113 114 115 116 117 118 119 120 121 122
      TYPE(t_potden),    INTENT(in)           :: vTot
      TYPE(t_oneD),      INTENT(in)           :: oneD
      TYPE(t_results),   INTENT(in)           :: results
      TYPE(t_potden),    INTENT(inout)        :: EnergyDen

      ! local
      INTEGER                         :: jspin

      TYPE(t_regionCharges)           :: regCharges
      TYPE(t_dos)                     :: dos
      TYPE(t_moments)                 :: moments
      TYPE(t_results)                 :: tmp_results
      TYPE(t_cdnvalJob)               :: cdnvalJob
      TYPE(t_potden)                  :: aux_den, real_den
123

Matthias Redies's avatar
Matthias Redies committed
124 125

      CALL regCharges%init(input, atoms)
126
      CALL dos%init(input,        atoms, kpts, vacuum)
127 128
!      CALL moments%init(input,    atoms)
      CALL moments%init(mpi,input,sphhar,atoms)
Matthias Redies's avatar
Matthias Redies committed
129 130 131 132
      tmp_results = results

      DO jspin = 1,input%jspins
         CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin)
133

Matthias Redies's avatar
Matthias Redies committed
134 135 136 137 138

         ! replace brillouin weights with auxillary weights
         CALL calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, cdnvalJob%weights)

         CALL cdnval(eig_id, mpi, kpts, jspin, noco, input, banddos, cell, atoms, &
139
            enpara, stars, vacuum,  sphhar, sym, vTot, oneD, cdnvalJob, &
140
            EnergyDen, regCharges, dos, tmp_results, moments, gfinp, hub1inp)
Matthias Redies's avatar
Matthias Redies committed
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
      ENDDO

   END SUBROUTINE calc_EnergyDen

   SUBROUTINE calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, f_ik)
      USE m_types_kpts
      USE m_eig66_io
      IMPLICIT NONE
      ! calculates new (auxillary-)weights as
      ! f_iks = w_iks * E_iks
      !, where  f_iks are the new (auxillary-)weights
      ! w_iks are the weights used in brillouin zone integration
      ! E_iks are the eigen energies

      INTEGER,      INTENT(in)        :: eig_id
      INTEGER,      INTENT(in)        :: jspin
      TYPE(t_kpts), INTENT(in)        :: kpts
      REAL,         INTENT(inout)     :: f_ik(:,:) ! f_ik(band_idx, kpt_idx)

      ! local vars
161
      REAL                       :: e_i(SIZE(f_ik,dim=1))
Matthias Redies's avatar
Matthias Redies committed
162 163 164
      INTEGER                    :: ikpt

      DO ikpt = 1,kpts%nkpt
165 166
         CALL read_eig(eig_id,ikpt,jspin, eig=e_i)
         f_ik(:,ikpt) = f_ik(:,ikpt) * e_i
Matthias Redies's avatar
Matthias Redies committed
167 168
      ENDDO
   END SUBROUTINE calc_EnergyDen_auxillary_weights
Matthias Redies's avatar
Matthias Redies committed
169

Matthias Redies's avatar
Matthias Redies committed
170
   subroutine set_zPrime(dim_idx, zMat, kpt, lapw, cell, zPrime)
Matthias Redies's avatar
Matthias Redies committed
171
      USE m_types
Matthias Redies's avatar
Matthias Redies committed
172
      USE m_constants
Matthias Redies's avatar
Matthias Redies committed
173
      implicit none
Matthias Redies's avatar
Matthias Redies committed
174
      INTEGER, intent(in)      :: dim_idx
Matthias Redies's avatar
Matthias Redies committed
175
      TYPE (t_mat), intent(in) :: zMat
176
      REAL, intent(in)         :: kpt(3)
Matthias Redies's avatar
Matthias Redies committed
177
      TYPE(t_lapw), intent(in) :: lapw
Matthias Redies's avatar
Matthias Redies committed
178
      TYPE(t_cell), intent(in) :: cell
Matthias Redies's avatar
Matthias Redies committed
179 180
      TYPE (t_mat)             :: zPrime

Matthias Redies's avatar
Matthias Redies committed
181
      REAL                     :: k_plus_g(3), fac
Matthias Redies's avatar
Matthias Redies committed
182
      INTEGER                  :: basis_idx
Matthias Redies's avatar
Matthias Redies committed
183

Matthias Redies's avatar
Matthias Redies committed
184 185
      call zPrime%free()
      call zPrime%init(zMat)
Matthias Redies's avatar
Matthias Redies committed
186

Matthias Redies's avatar
Matthias Redies committed
187
      do basis_idx = 1,size(lapw%gvec,dim=2)
Matthias Redies's avatar
Matthias Redies committed
188 189 190 191
         k_plus_g = kpt + lapw%gvec(:,basis_idx,1)
         k_plus_g = internal_to_rez(cell, k_plus_g)

         fac = k_plus_g(dim_idx)
Matthias Redies's avatar
Matthias Redies committed
192
         if(zPrime%l_real) then
193
            zPrime%data_r(basis_idx,:) =            fac * zMat%data_r(basis_idx,:)
Matthias Redies's avatar
Matthias Redies committed
194
         else
195
            zPrime%data_c(basis_idx,:) = ImagUnit * fac * zMat%data_c(basis_idx,:)
Matthias Redies's avatar
Matthias Redies committed
196
         endif
Matthias Redies's avatar
Matthias Redies committed
197
      enddo
Matthias Redies's avatar
Matthias Redies committed
198
   end subroutine set_zPrime
Matthias Redies's avatar
Matthias Redies committed
199 200 201 202 203 204 205 206 207 208

   function internal_to_rez(cell, vec) result(res)
      use m_types
      implicit none
      type(t_cell), intent(in) :: cell
      real, intent(in)      :: vec(3)
      real                  :: res(3)

      res = matmul(transpose(cell%bmat), vec)
   end function internal_to_rez
209

210 211 212 213 214 215 216 217 218
   subroutine undo_vgen_finalize(vtot, atoms, noco, stars)
      use m_types
      use m_constants
      use m_judft
      implicit none
      TYPE(t_potden), intent(inout)  :: vtot
      type(t_atoms), intent(in)      :: atoms
      type(t_noco), intent(in)       :: noco
      type(t_stars), intent(in)      :: stars
219

220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
      integer                        :: js, n, st

      do js = 1,size(vtot%mt,4)
         do n = 1,atoms%ntype
            vTot%mt(:atoms%jri(n),0,n,js) = vtot%mt(:atoms%jri(n),0,n,js) &
                  / (atoms%rmsh(:atoms%jri(n),n) / sfp_const )
         enddo
      enddo

      if(.not. noco%l_noco) then
         do js=1,size(vtot%pw_w,2)
            do st=1,stars%ng3
               vTot%pw_w(st,js) = vTot%pw_w(st,js) * stars%nstr(st)
            enddo
         enddo
      else
         call juDFT_error("undo vgen_finalize not implemented for noco")
      endif
   end subroutine undo_vgen_finalize

240
   subroutine set_kinED(mpi,   sphhar, atoms, sym,  xcpot, &
241
                        input, noco,   stars, vacuum,oned,cell,     den,     EnergyDen, vTot,kinED)
242
      use m_types
243
      use m_cdn_io
244 245 246 247
      implicit none
      TYPE(t_mpi),INTENT(IN)       :: mpi
      TYPE(t_sphhar),INTENT(IN)    :: sphhar
      TYPE(t_atoms),INTENT(IN)     :: atoms
248
      TYPE(t_sym), INTENT(IN)      :: sym
249
      CLASS(t_xcpot),INTENT(IN)    :: xcpot
250 251 252
      TYPE(t_input),INTENT(IN)     :: input
      TYPE(t_noco),INTENT(IN)      :: noco
      TYPE(t_stars),INTENT(IN)     :: stars
253 254
      TYPE(t_vacuum),INTENT(IN)    :: vacuum
      TYPE(t_oneD),INTENT(IN)      :: oneD
255 256
      TYPE(t_cell),INTENT(IN)      :: cell
      TYPE(t_potden),INTENT(IN)    :: den, EnergyDen, vTot
257
      TYPE(t_kinED),INTENT(OUT)    :: kinED
258

259
      TYPE(t_potden)               :: vTot_corrected
260 261 262 263 264 265 266 267

      LOGICAL :: perform_MetaGGA
      TYPE(t_potden)    :: core_den, val_den
      real :: rdum
      logical :: ldum
      perform_MetaGGA = ALLOCATED(EnergyDen%mt) &
                      .AND. (xcpot%exc_is_MetaGGA() .or. xcpot%vx_is_MetaGGA())
      if(.not.perform_MetaGGA) return
268
#ifdef CPP_LIBXC
269 270 271

      call readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,&
                          CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
272
                           0,rdum,ldum,core_den,'cdnc')
273 274
      call val_den%subPotDen(den,core_den)

275 276
      call vTot_corrected%copyPotDen(vTot)
      call undo_vgen_finalize(vTot_corrected, atoms, noco, stars)
277

278
      call set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot_corrected,kinED)
279
      call set_kinED_mt(mpi,   sphhar,    atoms, sym, noco,core_den, val_den, &
280
                           xcpot, EnergyDen, input, vTot_corrected,kinED)
281
#endif
282
   end subroutine set_kinED
283
#ifdef CPP_LIBXC
284
   subroutine set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot,kinED)
285 286 287
      use m_types
      use m_pw_tofrom_grid
      implicit none
288
      CLASS(t_xcpot),INTENT(IN)    :: xcpot
289 290 291
      TYPE(t_input),INTENT(IN)     :: input
      TYPE(t_noco),INTENT(IN)      :: noco
      TYPE(t_stars),INTENT(IN)     :: stars
292
      TYPE(t_sym), INTENT(IN)      :: sym
293 294
      TYPE(t_cell),INTENT(IN)      :: cell
      TYPE(t_potden),INTENT(IN)    :: den, EnergyDen, vTot
295
      TYPE(t_kinED),INTENT(INOUT)  :: kinED
296 297 298 299

      !local arrays
      REAL, ALLOCATABLE            :: den_rs(:,:), ED_rs(:,:), vTot_rs(:,:)
      TYPE(t_gradients)            :: tmp_grad
300

301
      CALL init_pw_grid(xcpot%needs_grad(),stars,sym,cell)
302

303 304 305 306 307 308
      CALL pw_to_grid(xcpot%needs_grad(), input%jspins, noco%l_noco, stars, &
                      cell,  EnergyDen%pw, tmp_grad, xcpot,    ED_rs)
      CALL pw_to_grid(xcpot%needs_grad(), input%jspins, noco%l_noco, stars, &
                      cell,  vTot%pw,      tmp_grad, xcpot,    vTot_rs)
      CALL pw_to_grid(xcpot%needs_grad(), input%jspins, noco%l_noco, stars, &
                      cell,  den%pw,       tmp_grad, xcpot,   den_rs)
309

310
      CALL finish_pw_grid()
311

312
      call calc_kinEnergyDen_pw(ED_rs, vTot_rs, den_rs, kinED%is)
Matthias Redies's avatar
Matthias Redies committed
313
      !xcpot%kinED%is  = ED_RS - vTot_RS * den_RS
314
      kinED%set = .True.
315 316
   end subroutine set_kinED_is

317
   subroutine set_kinED_mt(mpi,   sphhar,    atoms, sym, noco,core_den, val_den, &
318
                           xcpot, EnergyDen, input, vTot,kinED)
319 320 321 322 323 324
      use m_types
      use m_mt_tofrom_grid
      implicit none
      TYPE(t_mpi),INTENT(IN)         :: mpi
      TYPE(t_sphhar),INTENT(IN)      :: sphhar
      TYPE(t_atoms),INTENT(IN)       :: atoms
325
      TYPE(t_sym), INTENT(IN)        :: sym
326
      TYPE(t_noco), INTENT(IN)       :: noco
327
      TYPE(t_potden),INTENT(IN)      :: core_den, val_den, EnergyDen, vTot
328
      CLASS(t_xcpot),INTENT(IN)      :: xcpot
329
      TYPE(t_input),INTENT(IN)       :: input
330
      TYPE(t_kinED),INTENT(INOUT)    :: kinED
331 332 333 334 335 336 337 338 339 340 341 342 343
      INTEGER                        :: jr, loc_n, n, n_start, n_stride, cnt
      REAL,ALLOCATABLE               :: vTot_mt(:,:,:), ED_rs(:,:), vTot_rs(:,:), vTot0_rs(:,:),&
                                        core_den_rs(:,:), val_den_rs(:,:)
      TYPE(t_gradients)              :: tmp_grad
      TYPE(t_sphhar)                 :: tmp_sphhar

#ifdef CPP_MPI
      n_start=mpi%irank+1
      n_stride=mpi%isize
#else
      n_start=1
      n_stride=1
#endif
344
      CALL init_mt_grid(input%jspins,atoms,sphhar,xcpot%needs_grad(),sym)
345
      loc_n = 0
346 347 348 349 350 351
      allocate(ED_rs(atoms%nsp()*atoms%jmtd, input%jspins))
      allocate(vTot_rs, mold=ED_rs)
      allocate(vTot0_rs, mold=ED_rs)
      allocate(core_den_rs, mold=ED_rs)
      allocate(val_den_rs, mold=ED_rs)

352
      call kinED%alloc_mt(atoms%nsp()*atoms%jmtd, input%jspins, &
353 354 355 356 357 358 359 360 361 362
                                n_start,                atoms%ntype,  n_stride)
      loc_n = 0
      do n = n_start,atoms%ntype,n_stride
         loc_n = loc_n + 1

         if(.not. allocated(vTot_mt)) then
            allocate(vTot_mt(lbound(vTot%mt, dim=1):ubound(vTot%mt, dim=1),&
                             lbound(vTot%mt, dim=2):ubound(vTot%mt, dim=2),&
                             lbound(vTot%mt, dim=4):ubound(vTot%mt, dim=4)))
         endif
363

364 365 366
         do jr=1,atoms%jri(n)
            vTot_mt(jr,0:,:) = vTot%mt(jr,0:,n,:) * atoms%rmsh(jr,n)**2
         enddo
367
         CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms,sym, sphhar,.TRUE., EnergyDen%mt(:, 0:, n, :), &
368
                         n,  noco,   tmp_grad,     ED_rs)
369
         CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,sphhar,.TRUE., vTot_mt(:,0:,:), &
370
                         n,     noco,tmp_grad,     vTot_rs)
371

372 373 374
         tmp_sphhar%nlhd = sphhar%nlhd
         tmp_sphhar%nlh  = [(0, cnt=1,size(sphhar%nlh))]

375
         CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,tmp_sphhar,.TRUE., vTot_mt(:,0:0,:), &
376
                         n,    noco, tmp_grad,     vTot0_rs)
377
         CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,sphhar,.TRUE., &
378
                         core_den%mt(:,0:,n,:), n,noco, tmp_grad, core_den_rs)
379
         CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,sphhar,.TRUE., &
380
                         val_den%mt(:,0:,n,:), n,noco, tmp_grad, val_den_rs)
381

Matthias Redies's avatar
Matthias Redies committed
382
         call calc_kinEnergyDen_mt(ED_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, &
383
                                   kinED%mt(:,:,loc_n))
Matthias Redies's avatar
Matthias Redies committed
384
         !xcpot%kinED%mt(:,:,loc_n) = ED_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
385
      enddo
386
      kinED%set = .True.
387
      CALL finish_mt_grid()
388
   end subroutine set_kinED_mt
389
#endif
Matthias Redies's avatar
Matthias Redies committed
390
END MODULE m_metagga