metagga.F90 8.69 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
    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 13 14

    type t_RS_potden
       REAL, ALLOCATABLE :: is(:,:), mt(:,:)
    end type t_RS_potden

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

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

Matthias Redies's avatar
Matthias Redies committed
30
        !implicit allocation
31
        kinEnergyDen_RS = EnergyDen_RS - vTot_RS * den_RS
Matthias Redies's avatar
cut off  
Matthias Redies committed
32 33 34 35 36 37 38

        if(any(kinEnergyDen_RS < eps)) then
           write (6,*) "         lowest kinetic energy density cutoff = ", minval(kinEnergyDen_RS)
           kinEnergyDen_RS = max(kinEnergyDen_RS, eps)
           write (6,*) "         afterwards = ", minval(kinEnergyDen_RS)
        endif

Matthias Redies's avatar
Matthias Redies committed
39
#else
40 41
        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
42
#endif
43
    END SUBROUTINE calc_kinEnergyDen
Matthias Redies's avatar
Matthias Redies committed
44

45 46
    SUBROUTINE calc_EnergyDen(eig_id, mpi, kpts, noco, input, banddos, cell, atoms, enpara, stars, &
                              vacuum, DIMENSION, sphhar, sym, vTot, oneD, results, EnergyDen)
Matthias Redies's avatar
Matthias Redies committed
47 48 49 50 51 52
        ! 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
         
        
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
        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
Matthias Redies's avatar
Matthias Redies committed
84 85 86


        ! local
87
        INTEGER                         :: jspin
Matthias Redies's avatar
Matthias Redies committed
88
 
89 90 91 92 93 94 95 96 97 98
        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)
Matthias Redies's avatar
Matthias Redies committed
99 100
        tmp_results = results

101 102
        DO jspin = 1,input%jspins
            CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin)
Matthias Redies's avatar
Matthias Redies committed
103 104

            ! replace brillouin weights with auxillary weights
105
            CALL calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, cdnvalJob%weights)
Matthias Redies's avatar
Matthias Redies committed
106

107 108
            CALL cdnval(eig_id, mpi, kpts, jspin, noco, input, banddos, cell, atoms, &
                        enpara, stars, vacuum, DIMENSION, sphhar, sym, vTot, oneD, cdnvalJob, &
109
                        EnergyDen, regCharges, dos, tmp_results, moments)
110
        ENDDO
Matthias Redies's avatar
Matthias Redies committed
111

112
    END SUBROUTINE calc_EnergyDen
Matthias Redies's avatar
Matthias Redies committed
113

114 115 116 117
    SUBROUTINE calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, f_ik)
        USE m_types_kpts
        USE m_eig66_io
        IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
118 119 120 121 122 123 124
        ! 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


125 126 127 128
        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)
Matthias Redies's avatar
Matthias Redies committed
129 130

        ! local vars
131 132
        REAL                       :: w_i(SIZE(f_ik,dim=1)), e_i(SIZE(f_ik,dim=1))
        INTEGER                    :: ikpt
Matthias Redies's avatar
Matthias Redies committed
133

134 135
        DO ikpt = 1,kpts%nkpt
            CALL read_eig(eig_id,ikpt,jspin, eig=e_i, w_iks=w_i)
Matthias Redies's avatar
Matthias Redies committed
136 137

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

Matthias Redies's avatar
Matthias Redies committed
141 142 143 144 145 146 147 148
    !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
149

Matthias Redies's avatar
Matthias Redies committed
150
        !IMPLICIT NONE 
Matthias Redies's avatar
Matthias Redies committed
151
        
Matthias Redies's avatar
Matthias Redies committed
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
        !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
Matthias Redies's avatar
Matthias Redies committed
199

200
    FUNCTION subtract_RS(rs1, rs2) RESULT(rs_out)
Matthias Redies's avatar
Matthias Redies committed
201
       USE m_types_xcpot
202
        IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
203

Matthias Redies's avatar
Matthias Redies committed
204 205
        TYPE(t_RS_potden), INTENT(in)  :: rs1, rs2
        TYPE(t_RS_potden)              :: rs_out
Matthias Redies's avatar
Matthias Redies committed
206 207

        rs_out%is = rs1%is - rs2%is 
Matthias Redies's avatar
Matthias Redies committed
208
        rs_out%mt = rs1%mt - rs2%mt
209
    END FUNCTION subtract_RS 
Matthias Redies's avatar
Matthias Redies committed
210

211
    FUNCTION multiply_RS(rs1, rs2) RESULT(rs_out)
Matthias Redies's avatar
Matthias Redies committed
212
       USE m_types_xcpot
213
        IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
214
        
Matthias Redies's avatar
Matthias Redies committed
215 216
        TYPE(t_RS_potden), INTENT(in)  :: rs1, rs2
        TYPE(t_RS_potden)              :: rs_out
Matthias Redies's avatar
Matthias Redies committed
217

Matthias Redies's avatar
Matthias Redies committed
218 219
        rs_out%is = rs1%is * rs2%is
        rs_out%mt = rs1%mt * rs2%mt
220
    END FUNCTION multiply_RS
Matthias Redies's avatar
Matthias Redies committed
221

222
END MODULE m_metagga