metagga.F90 8.04 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 8
   PUBLIC  :: calc_EnergyDen
   PRIVATE :: calc_EnergyDen_auxillary_weights, subtract_RS, multiply_RS
Matthias Redies's avatar
Matthias Redies committed
9

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

Matthias Redies's avatar
Matthias Redies committed
14 15 16
   INTERFACE OPERATOR (-)
      PROCEDURE subtract_RS
   END INTERFACE OPERATOR (-)
17

Matthias Redies's avatar
Matthias Redies committed
18 19 20
   INTERFACE OPERATOR (*)
      PROCEDURE multiply_RS
   END INTERFACE OPERATOR (*)
21
CONTAINS
Matthias Redies's avatar
Matthias Redies committed
22 23 24 25 26
   SUBROUTINE calc_kinEnergyDen(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS)
#ifdef CPP_LIBXC
      IMPLICIT NONE
      REAL, INTENT(in)                 :: den_RS(:,:), EnergyDen_RS(:,:), vTot_RS(:,:)
      REAL, INTENT(inout), allocatable :: kinEnergyDen_RS(:,:)
Matthias Redies's avatar
Matthias Redies committed
27
      REAL, PARAMETER                  :: eps = 1e-15
Matthias Redies's avatar
Matthias Redies committed
28 29 30 31 32 33 34 35

      !implicit allocation
      kinEnergyDen_RS = EnergyDen_RS - vTot_RS * den_RS

      if(any(kinEnergyDen_RS < eps)) then
         write (6,*) "         lowest kinetic energy density cutoff = ", minval(kinEnergyDen_RS)
         kinEnergyDen_RS = max(kinEnergyDen_RS, eps)
      endif
Matthias Redies's avatar
Matthias Redies committed
36
   
Matthias Redies's avatar
Matthias Redies committed
37
#else
Matthias Redies's avatar
Matthias Redies committed
38 39
      USE m_juDFT_stop
      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
40
#endif
Matthias Redies's avatar
Matthias Redies committed
41 42 43 44 45 46 47 48 49 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 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 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 212 213
   END SUBROUTINE calc_kinEnergyDen

   SUBROUTINE calc_EnergyDen(eig_id, mpi, kpts, noco, input, banddos, cell, atoms, enpara, stars, &
                             vacuum, DIMENSION, sphhar, sym, vTot, oneD, results, EnergyDen)
      ! 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

      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)

         ! 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, &
                     enpara, stars, vacuum, DIMENSION, sphhar, sym, vTot, oneD, cdnvalJob, &
                     EnergyDen, regCharges, dos, tmp_results, moments)
      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
      REAL                       :: w_i(SIZE(f_ik,dim=1)), e_i(SIZE(f_ik,dim=1))
      INTEGER                    :: ikpt

      DO ikpt = 1,kpts%nkpt
         CALL read_eig(eig_id,ikpt,jspin, eig=e_i, w_iks=w_i)

         f_ik(:,ikpt) = e_i * w_i
      ENDDO
   END SUBROUTINE calc_EnergyDen_auxillary_weights

   !SUBROUTINE transform_to_grid(input, noco, sym, stars, cell, den, atoms, sphhar, EnergyDen, vTot, den_RS, EnergyDen_RS, vTot_RS)
   !USE m_types_potden
   !USE m_types_setup
   !USE m_types_xcpot_libxc
   !USE m_types_xcpot
   !USE m_juDFT_stop
   !USE m_pw_tofrom_grid
   !USE m_mt_tofrom_grid

   !IMPLICIT NONE

   !TYPE(t_potden),    INTENT(in)        :: den, EnergyDen, vTot
   !TYPE(t_input),     INTENT(in)        :: input
   !TYPE(t_noco),      INTENT(in)        :: noco
   !TYPE(t_sym),       INTENT(in)        :: sym
   !TYPE(t_stars),     INTENT(in)        :: stars
   !TYPE(t_cell),      INTENT(in)        :: cell
   !TYPE(t_atoms),     INTENT(in)        :: atoms
   !TYPE(t_sphhar),    INTENT(in)        :: sphhar
   !REAL(:,:), INTENT(out)               :: den_RS, EnergyDen_RS, vTot_RS ! could be changed to a real-space type

   !!local vars
   !TYPE(t_xcpot_libxc) ::aux_xcpot
   !TYPE(t_gradients)   :: tmp_grad
   !INTEGER, PARAMETER  :: id_corr = 9, id_exch = 1
   !INTEGER             :: nsp, n

   !!make some auxillary xcpot, that is not a GGA (we don't need gradients)
   !CALL aux_xcpot%init(input%jspins, id_exch, id_corr, id_exch, id_corr)
   !IF(aux_xcpot%vxc_is_gga()) CALL juDFT_error("aux_xcpot must not be GGA", &
   !hint="choose id_corr and id_exch correctly")

   !! interstitial part
   !CALL init_pw_grid(aux_xcpot,stars,sym,cell)

   !CALL pw_to_grid(aux_xcpot, input%jspins, noco%l_noco, stars, cell, den%pw,       tmp_grad, den_RS%is)
   !CALL pw_to_grid(aux_xcpot, input%jspins, noco%l_noco, stars, cell, EnergyDen%pw, tmp_grad, EnergyDen_RS%is)
   !CALL pw_to_grid(aux_xcpot, input%jspins, noco%l_noco, stars, cell, vTot%pw,      tmp_grad, vTot_RS%is)

   !CALL finish_pw_grid()

   !! muffin tins
   !nsp=(atoms%lmaxd+1+MOD(atoms%lmaxd+1,2))*(2*atoms%lmaxd+1)
   !CALL init_mt_grid(nsp,input%jspins,atoms,sphhar,aux_xcpot,sym)

   !DO n = 1,atoms%ntype
   !CALL mt_to_grid(aux_xcpot, input%jspins, atoms,    sphhar, den%mt(:,0:,n,:), &
   !nsp,       n,            tmp_grad, den_RS%mt)
   !CALL mt_to_grid(aux_xcpot, input%jspins, atoms,    sphhar, EnergyDen%mt(:,0:,n,:), &
   !nsp,       n,            tmp_grad, EnergyDen_RS%mt)
   !CALL mt_to_grid(aux_xcpot, input%jspins, atoms,    sphhar, vTot%mt(:,0:,n,:), &
   !nsp,       n,            tmp_grad, vTot_RS%mt)
   !ENDDO

   !CALL finish_mt_grid()
   !END SUBROUTINE transform_to_grid

   FUNCTION subtract_RS(rs1, rs2) RESULT(rs_out)
      USE m_types_xcpot
      IMPLICIT NONE

      TYPE(t_RS_potden), INTENT(in)  :: rs1, rs2
      TYPE(t_RS_potden)              :: rs_out

      rs_out%is = rs1%is - rs2%is
      rs_out%mt = rs1%mt - rs2%mt
   END FUNCTION subtract_RS

   FUNCTION multiply_RS(rs1, rs2) RESULT(rs_out)
      USE m_types_xcpot
      IMPLICIT NONE

      TYPE(t_RS_potden), INTENT(in)  :: rs1, rs2
      TYPE(t_RS_potden)              :: rs_out

      rs_out%is = rs1%is * rs2%is
      rs_out%mt = rs1%mt * rs2%mt
   END FUNCTION multiply_RS
Matthias Redies's avatar
Matthias Redies committed
214

215
END MODULE m_metagga