metagga.F90 13 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
CONTAINS
17
   SUBROUTINE calc_kinEnergyDen_pw(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS)
18
      USE m_juDFT_stop
Matthias Redies's avatar
Matthias Redies committed
19
      !use m_cdngen
Matthias Redies's avatar
Matthias Redies committed
20 21 22
      IMPLICIT NONE
      REAL, INTENT(in)                 :: den_RS(:,:), EnergyDen_RS(:,:), vTot_RS(:,:)
      REAL, INTENT(inout), allocatable :: kinEnergyDen_RS(:,:)
23
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
24
      REAL, PARAMETER                  :: eps = 1e-15
Matthias Redies's avatar
Matthias Redies committed
25 26

      kinEnergyDen_RS = EnergyDen_RS - vTot_RS * den_RS
27 28 29 30 31 32
#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
33
                                   kinEnergyDen_RS)
34 35 36 37
      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(:,:)
38
      REAL, INTENT(inout)              :: kinEnergyDen_RS(:,:)
39 40 41

#ifdef CPP_LIBXC
      kinEnergyDen_RS = EnergyDen_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
Matthias Redies's avatar
Matthias Redies committed
42
#else
Matthias Redies's avatar
Matthias Redies committed
43
      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
44
#endif
45
   END SUBROUTINE calc_kinEnergyDen_mt
Matthias Redies's avatar
Matthias Redies committed
46

Matthias Redies's avatar
Matthias Redies committed
47

Matthias Redies's avatar
Matthias Redies committed
48
   SUBROUTINE calc_EnergyDen(eig_id, mpi, kpts, noco, input, banddos, cell, atoms, enpara, stars, &
49
         vacuum, DIMENSION, sphhar, sym, vTot, oneD, results, EnergyDen)
Matthias Redies's avatar
Matthias Redies committed
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
      ! 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
      TYPE(t_dimension), INTENT(in)           :: DIMENSION
      TYPE(t_sphhar),    INTENT(in)           :: sphhar
      TYPE(t_sym),       INTENT(in)           :: sym
      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
96

Matthias Redies's avatar
Matthias Redies committed
97 98 99 100 101 102 103 104

      CALL regCharges%init(input, atoms)
      CALL dos%init(input,        atoms, DIMENSION, kpts, vacuum)
      CALL moments%init(input,    atoms)
      tmp_results = results

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

Matthias Redies's avatar
Matthias Redies committed
106 107 108 109 110

         ! 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, &
111 112
            enpara, stars, vacuum, DIMENSION, sphhar, sym, vTot, oneD, cdnvalJob, &
            EnergyDen, regCharges, dos, tmp_results, moments)
Matthias Redies's avatar
Matthias Redies committed
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
      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
133
      REAL                       :: e_i(SIZE(f_ik,dim=1))
Matthias Redies's avatar
Matthias Redies committed
134 135 136
      INTEGER                    :: ikpt

      DO ikpt = 1,kpts%nkpt
137 138
         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
139 140
      ENDDO
   END SUBROUTINE calc_EnergyDen_auxillary_weights
Matthias Redies's avatar
Matthias Redies committed
141

Matthias Redies's avatar
Matthias Redies committed
142
   subroutine set_zPrime(dim_idx, zMat, kpt, lapw, cell, zPrime)
Matthias Redies's avatar
Matthias Redies committed
143
      USE m_types
Matthias Redies's avatar
Matthias Redies committed
144
      USE m_constants
Matthias Redies's avatar
Matthias Redies committed
145
      implicit none
Matthias Redies's avatar
Matthias Redies committed
146
      INTEGER, intent(in)      :: dim_idx
Matthias Redies's avatar
Matthias Redies committed
147 148 149
      TYPE (t_mat), intent(in) :: zMat
      REAL, intent(in)         :: kpt(3) 
      TYPE(t_lapw), intent(in) :: lapw
Matthias Redies's avatar
Matthias Redies committed
150
      TYPE(t_cell), intent(in) :: cell
Matthias Redies's avatar
Matthias Redies committed
151 152
      TYPE (t_mat)             :: zPrime

Matthias Redies's avatar
Matthias Redies committed
153
      REAL                     :: k_plus_g(3), fac
Matthias Redies's avatar
Matthias Redies committed
154
      INTEGER                  :: basis_idx
Matthias Redies's avatar
Matthias Redies committed
155

Matthias Redies's avatar
Matthias Redies committed
156 157
      call zPrime%free()
      call zPrime%init(zMat)
Matthias Redies's avatar
Matthias Redies committed
158

Matthias Redies's avatar
Matthias Redies committed
159
      do basis_idx = 1,size(lapw%gvec,dim=2)
Matthias Redies's avatar
Matthias Redies committed
160 161 162 163
         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
164
         if(zPrime%l_real) then
Matthias Redies's avatar
Matthias Redies committed
165
            zPrime%data_r(basis_idx,:) =            fac * zMat%data_r(basis_idx,:) 
Matthias Redies's avatar
Matthias Redies committed
166
         else
Matthias Redies's avatar
Matthias Redies committed
167
            zPrime%data_c(basis_idx,:) = ImagUnit * fac * zMat%data_c(basis_idx,:) 
Matthias Redies's avatar
Matthias Redies committed
168
         endif
Matthias Redies's avatar
Matthias Redies committed
169
      enddo
Matthias Redies's avatar
Matthias Redies committed
170
   end subroutine set_zPrime
Matthias Redies's avatar
Matthias Redies committed
171 172 173 174 175 176 177 178 179 180

   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
181

182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
   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
   
      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

212
   subroutine set_kinED(mpi,   sphhar, atoms, sym, core_den, val_den, xcpot, &
213 214 215 216 217 218
                        input, noco,   stars, cell,     den,     EnergyDen, vTot)
      use m_types
      implicit none
      TYPE(t_mpi),INTENT(IN)       :: mpi
      TYPE(t_sphhar),INTENT(IN)    :: sphhar
      TYPE(t_atoms),INTENT(IN)     :: atoms
219
      TYPE(t_sym), INTENT(IN)      :: sym
220 221 222 223 224 225 226
      TYPE(t_potden),INTENT(IN)    :: core_den, val_den
      CLASS(t_xcpot),INTENT(INOUT) :: xcpot
      TYPE(t_input),INTENT(IN)     :: input
      TYPE(t_noco),INTENT(IN)      :: noco
      TYPE(t_stars),INTENT(IN)     :: stars
      TYPE(t_cell),INTENT(IN)      :: cell
      TYPE(t_potden),INTENT(IN)    :: den, EnergyDen, vTot
227 228 229 230 231
      
      TYPE(t_potden)               :: vTot_corrected
   
      call vTot_corrected%copyPotDen(vTot)
      call undo_vgen_finalize(vTot_corrected, atoms, noco, stars)
232

233
      call set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot_corrected)
234
      call set_kinED_mt(mpi,   sphhar,    atoms, sym, core_den, val_den, &
235
                           xcpot, EnergyDen, input, vTot_corrected)
236 237
   end subroutine set_kinED

238
   subroutine set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot)
239 240 241 242 243 244 245
      use m_types
      use m_pw_tofrom_grid
      implicit none
      CLASS(t_xcpot),INTENT(INOUT) :: xcpot
      TYPE(t_input),INTENT(IN)     :: input
      TYPE(t_noco),INTENT(IN)      :: noco
      TYPE(t_stars),INTENT(IN)     :: stars
246
      TYPE(t_sym), INTENT(IN)      :: sym
247 248 249 250 251 252
      TYPE(t_cell),INTENT(IN)      :: cell
      TYPE(t_potden),INTENT(IN)    :: den, EnergyDen, vTot

      !local arrays
      REAL, ALLOCATABLE            :: den_rs(:,:), ED_rs(:,:), vTot_rs(:,:)
      TYPE(t_gradients)            :: tmp_grad
253 254
      
      CALL init_pw_grid(xcpot,stars,sym,cell)
255 256 257 258 259 260 261 262

      CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
                      cell,  EnergyDen%pw, tmp_grad,    ED_rs)
      CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
                      cell,  vTot%pw,      tmp_grad,    vTot_rs)
      CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
                      cell,  den%pw,       tmp_grad,    den_rs)

263
      CALL finish_pw_grid()
Matthias Redies's avatar
Matthias Redies committed
264 265 266
      
      call calc_kinEnergyDen_pw(ED_rs, vTot_rs, den_rs, xcpot%kinED%is)
      !xcpot%kinED%is  = ED_RS - vTot_RS * den_RS
267 268 269
      xcpot%kinED%set = .True.
   end subroutine set_kinED_is

270
   subroutine set_kinED_mt(mpi,   sphhar,    atoms, sym, core_den, val_den, &
271 272 273 274 275 276 277
                           xcpot, EnergyDen, input, vTot)
      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
278
      TYPE(t_sym), INTENT(IN)        :: sym
279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295
      TYPE(t_potden),INTENT(IN)      :: core_den, val_den, EnergyDen, vTot
      CLASS(t_xcpot),INTENT(INOUT)   :: xcpot
      TYPE(t_input),INTENT(IN)       :: input

      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
296
      CALL init_mt_grid(input%jspins,atoms,sphhar,xcpot,sym)
297
      loc_n = 0
298 299 300 301 302 303
      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)

304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333
      call xcpot%kinED%alloc_mt(atoms%nsp()*atoms%jmtd, input%jspins, &
                                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
         
         do jr=1,atoms%jri(n)
            vTot_mt(jr,0:,:) = vTot%mt(jr,0:,n,:) * atoms%rmsh(jr,n)**2
         enddo
         CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, EnergyDen%mt(:, 0:, n, :), &
                         n,     tmp_grad,     ED_rs)
         CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, vTot_mt(:,0:,:), &
                         n,     tmp_grad,     vTot_rs)
         
         tmp_sphhar%nlhd = sphhar%nlhd
         tmp_sphhar%nlh  = [(0, cnt=1,size(sphhar%nlh))]

         CALL mt_to_grid(xcpot, input%jspins, atoms, tmp_sphhar, vTot_mt(:,0:0,:), &
                         n,     tmp_grad,     vTot0_rs)
         CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, &
                         core_den%mt(:,0:,n,:), n, tmp_grad, core_den_rs)
         CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, &
                         val_den%mt(:,0:,n,:), n, tmp_grad, val_den_rs)
         
Matthias Redies's avatar
Matthias Redies committed
334 335 336
         call calc_kinEnergyDen_mt(ED_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, &
                                   xcpot%kinED%mt(:,:,loc_n))
         !xcpot%kinED%mt(:,:,loc_n) = ED_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
337 338
      enddo
      xcpot%kinED%set = .True.
339
      CALL finish_mt_grid()
340
   end subroutine set_kinED_mt
Matthias Redies's avatar
Matthias Redies committed
341
END MODULE m_metagga