LCOV - code coverage report
Current view: top level - xc-pot - metagga.F90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 5.0 % 119 6
Test Date: 2025-06-14 04:34:23 Functions: 6.7 % 15 1

            Line data    Source code
       1              : !--------------------------------------------------------------------------------
       2              : ! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3              : ! This file is part of FLEUR and available as free software under the conditions
       4              : ! of the MIT license as expressed in the LICENSE file in more detail.
       5              : !--------------------------------------------------------------------------------
       6              : MODULE m_metagga
       7              :    PUBLIC  :: calc_EnergyDen
       8              :    PRIVATE :: calc_EnergyDen_auxillary_weights, &
       9              :               calc_kinEnergyDen_pw, &
      10              :               calc_kinEnergyDen_mt
      11              : 
      12              :    type t_RS_potden
      13              :       REAL, ALLOCATABLE :: is(:,:), mt(:,:)
      14              :    end type t_RS_potden
      15              : 
      16              :    TYPE t_kinED
      17              :       logical             :: set=.false.
      18              :       real, allocatable   :: is(:,:)   ! (nsp*jmtd, jspins)
      19              :       real, allocatable   :: mt(:,:,:) ! (nsp*jmtd, jspins, local num of types)
      20              :    contains
      21              :       procedure           :: alloc_mt => kED_alloc_mt
      22              :    END TYPE t_kinED
      23              : 
      24              : CONTAINS
      25              : 
      26            0 :   SUBROUTINE kED_alloc_mt(kED,nsp_x_jmtd, jspins, n_start, n_types, n_stride)
      27              :     IMPLICIT NONE
      28              :     class(t_kinED), intent(inout)   :: kED
      29              :     integer, intent(in)            :: nsp_x_jmtd, jspins, n_start, n_types, n_stride
      30              :     integer                        :: cnt, n
      31              : 
      32            0 :     if(.not. allocated(kED%mt)) then
      33            0 :       cnt = 0
      34            0 :       do n = n_start,n_types,n_stride
      35            0 :         cnt = cnt + 1
      36              :       enddo
      37            0 :       allocate(kED%mt(nsp_x_jmtd, jspins, cnt), source=0.0)
      38              :     endif
      39            0 :   END SUBROUTINE kED_alloc_mt
      40              : 
      41              : 
      42            0 :    SUBROUTINE calc_kinEnergyDen_pw(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS)
      43              :       USE m_juDFT_stop
      44              :       !use m_cdngen
      45              :       IMPLICIT NONE
      46              :       REAL, INTENT(in)                 :: den_RS(:,:), EnergyDen_RS(:,:), vTot_RS(:,:)
      47              :       REAL, INTENT(inout), allocatable :: kinEnergyDen_RS(:,:)
      48              : #ifdef CPP_LIBXC
      49              :       REAL, PARAMETER                  :: eps = 1e-15
      50              : 
      51            0 :       kinEnergyDen_RS = EnergyDen_RS - vTot_RS * den_RS
      52              : #else
      53              :       CALL juDFT_error("MetaGGA require LibXC",hint="compile Fleur with LibXC (e.g. by giving '-external libxc' to ./configure")
      54              : #endif
      55            0 :    END SUBROUTINE calc_kinEnergyDen_pw
      56              : 
      57            0 :    SUBROUTINE calc_kinEnergyDen_mt(EnergyDen_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, &
      58            0 :                                    kinEnergyDen_RS)
      59              :       USE m_juDFT_stop
      60              :       implicit none
      61              :       REAL, INTENT(in)                 :: EnergyDen_RS(:,:), vTot_rs(:,:), vTot0_rs(:,:), core_den_rs(:,:), val_den_rs(:,:)
      62              :       REAL, INTENT(inout)              :: kinEnergyDen_RS(:,:)
      63              : 
      64              : #ifdef CPP_LIBXC
      65            0 :       kinEnergyDen_RS = EnergyDen_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
      66              : #else
      67              :       CALL juDFT_error("MetaGGA require LibXC",hint="compile Fleur with LibXC (e.g. by giving '-external libxc' to ./configure")
      68              : #endif
      69            0 :    END SUBROUTINE calc_kinEnergyDen_mt
      70              : 
      71              : 
      72            0 :    SUBROUTINE calc_EnergyDen(eig_id, fmpi, kpts, noco,nococonv, input, banddos, cell, atoms, enpara, stars, &
      73              :          vacuum,  sphhar, sym, gfinp, hub1inp, vTot,   results, EnergyDen)
      74              :       ! calculates the energy density
      75              :       ! EnergyDen = \sum_i n_i(r) \varepsilon_i
      76              :       ! where n_i(r) is the one-particle density
      77              :       ! and \varepsilon_i are the eigenenergies
      78              : 
      79              :       USE m_types_setup
      80              :       USE m_types_potden
      81              :       USE m_types_kpts
      82              :       USE m_types_mpi
      83              :       USE m_types_enpara
      84              :       USE m_types_misc
      85              :       USE m_types_regionCharges
      86              :       USE m_types_dos
      87              :       USE m_types_vacdos
      88              :       USE m_types_cdnval
      89              :       USE m_cdnval
      90              :       use m_types_nococonv
      91              :       IMPLICIT NONE
      92              : 
      93              :       INTEGER,           INTENT(in)           :: eig_id
      94              :       TYPE(t_mpi),       INTENT(in)           :: fmpi
      95              :       TYPE(t_kpts),      INTENT(in)           :: kpts
      96              :       TYPE(t_noco),      INTENT(in)           :: noco
      97              :       TYPE(t_nococonv),  INTENT(in)           :: nococonv
      98              :       TYPE(t_input),     INTENT(in)           :: input
      99              :       TYPE(t_banddos),   INTENT(in)           :: banddos
     100              :       TYPE(t_cell),      INTENT(in)           :: cell
     101              :       TYPE(t_atoms),     INTENT(in)           :: atoms
     102              :       TYPE(t_enpara),    INTENT(in)           :: enpara
     103              :       TYPE(t_stars),     INTENT(in)           :: stars
     104              :       TYPE(t_vacuum),    INTENT(in)           :: vacuum
     105              : 
     106              :       TYPE(t_sphhar),    INTENT(in)           :: sphhar
     107              :       TYPE(t_sym),       INTENT(in)           :: sym
     108              :       TYPE(t_gfinp),     INTENT(in)           :: gfinp
     109              :       TYPE(t_hub1inp),   INTENT(in)           :: hub1inp
     110              :       TYPE(t_potden),    INTENT(in)           :: vTot
     111              :        
     112              :       TYPE(t_results),   INTENT(in)           :: results
     113              :       TYPE(t_potden),    INTENT(inout)        :: EnergyDen
     114              : 
     115              :       ! local
     116              :       INTEGER                         :: jspin
     117              : 
     118            0 :       TYPE(t_regionCharges)           :: regCharges
     119            0 :       TYPE(t_dos)                     :: dos
     120            0 :       TYPE(t_vacdos)                  :: vacdos
     121              : 
     122            0 :       TYPE(t_moments)                 :: moments
     123            0 :       TYPE(t_results)                 :: tmp_results
     124            0 :       TYPE(t_cdnvalJob)               :: cdnvalJob
     125              :       TYPE(t_potden)                  :: aux_den, real_den
     126              : 
     127              : 
     128            0 :       CALL regCharges%init(input, atoms)
     129            0 :       CALL dos%init(input,        atoms, kpts, banddos,results%eig)
     130            0 :       CALL vacdos%init(input,        atoms, kpts, banddos,results%eig)
     131              : !      CALL moments%init(input,    atoms)
     132            0 :       CALL moments%init(fmpi,input,sphhar,atoms)
     133            0 :       tmp_results = results
     134              : 
     135            0 :       DO jspin = 1,input%jspins
     136            0 :          CALL cdnvalJob%init(fmpi,input,kpts,noco,results,jspin)
     137              : 
     138              : 
     139              :          ! replace brillouin weights with auxillary weights
     140            0 :          CALL calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, cdnvalJob%weights)
     141              : 
     142              :          CALL cdnval(eig_id, fmpi, kpts, jspin, noco,nococonv, input, banddos, cell, atoms, &
     143              :             enpara, stars, vacuum,  sphhar, sym, vTot,   cdnvalJob, &
     144            0 :             EnergyDen, regCharges, dos, vacdos,tmp_results, moments, gfinp, hub1inp)
     145              :       ENDDO
     146              : 
     147            0 :    END SUBROUTINE calc_EnergyDen
     148              : 
     149            0 :    SUBROUTINE calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, f_ik)
     150              :       USE m_types_kpts
     151              :       USE m_eig66_io
     152              :       IMPLICIT NONE
     153              :       ! calculates new (auxillary-)weights as
     154              :       ! f_iks = w_iks * E_iks
     155              :       !, where  f_iks are the new (auxillary-)weights
     156              :       ! w_iks are the weights used in brillouin zone integration
     157              :       ! E_iks are the eigen energies
     158              : 
     159              :       INTEGER,      INTENT(in)        :: eig_id
     160              :       INTEGER,      INTENT(in)        :: jspin
     161              :       TYPE(t_kpts), INTENT(in)        :: kpts
     162              :       REAL,         INTENT(inout)     :: f_ik(:,:) ! f_ik(band_idx, kpt_idx)
     163              : 
     164              :       ! local vars
     165            0 :       REAL                       :: e_i(SIZE(f_ik,dim=1))
     166              :       INTEGER                    :: ikpt
     167              : 
     168            0 :       DO ikpt = 1,kpts%nkpt
     169            0 :          CALL read_eig(eig_id,ikpt,jspin, eig=e_i)
     170            0 :          f_ik(:,ikpt) = f_ik(:,ikpt) * e_i
     171              :       ENDDO
     172            0 :    END SUBROUTINE calc_EnergyDen_auxillary_weights
     173              : 
     174            0 :    subroutine set_zPrime(dim_idx, zMat, kpt, lapw, cell, zPrime)
     175              :       USE m_types
     176              :       USE m_constants
     177              :       implicit none
     178              :       INTEGER, intent(in)      :: dim_idx
     179              :       TYPE (t_mat), intent(in) :: zMat
     180              :       REAL, intent(in)         :: kpt(3)
     181              :       TYPE(t_lapw), intent(in) :: lapw
     182              :       TYPE(t_cell), intent(in) :: cell
     183              :       TYPE (t_mat)             :: zPrime
     184              : 
     185              :       REAL                     :: k_plus_g(3), fac
     186              :       INTEGER                  :: basis_idx
     187              : 
     188            0 :       call zPrime%free()
     189            0 :       call zPrime%init(zMat)
     190              : 
     191            0 :       do basis_idx = 1,size(lapw%gvec,dim=2)
     192            0 :          k_plus_g = kpt + lapw%gvec(:,basis_idx,1)
     193            0 :          k_plus_g = internal_to_rez(cell, k_plus_g)
     194              : 
     195            0 :          fac = k_plus_g(dim_idx)
     196            0 :          if(zPrime%l_real) then
     197            0 :             zPrime%data_r(basis_idx,:) =            fac * zMat%data_r(basis_idx,:)
     198              :          else
     199            0 :             zPrime%data_c(basis_idx,:) = ImagUnit * fac * zMat%data_c(basis_idx,:)
     200              :          endif
     201              :       enddo
     202            0 :    end subroutine set_zPrime
     203              : 
     204            0 :    function internal_to_rez(cell, vec) result(res)
     205              :       use m_types
     206              :       implicit none
     207              :       type(t_cell), intent(in) :: cell
     208              :       real, intent(in)      :: vec(3)
     209              :       real                  :: res(3)
     210              : 
     211            0 :       res = matmul(transpose(cell%bmat), vec)
     212            0 :    end function internal_to_rez
     213              : 
     214            0 :    subroutine undo_vgen_finalize(vtot, atoms, noco, stars)
     215              :       use m_types
     216              :       use m_constants
     217              :       use m_judft
     218              :       implicit none
     219              :       TYPE(t_potden), intent(inout)  :: vtot
     220              :       type(t_atoms), intent(in)      :: atoms
     221              :       type(t_noco), intent(in)       :: noco
     222              :       type(t_stars), intent(in)      :: stars
     223              : 
     224              :       integer                        :: js, n, st
     225              : 
     226            0 :       do js = 1,size(vtot%mt,4)
     227            0 :          do n = 1,atoms%ntype
     228              :             vTot%mt(:atoms%jri(n),0,n,js) = vtot%mt(:atoms%jri(n),0,n,js) &
     229            0 :                   / (atoms%rmsh(:atoms%jri(n),n) / sfp_const )
     230              :          enddo
     231              :       enddo
     232              : 
     233            0 :       if(.not. noco%l_noco) then
     234            0 :          do js=1,size(vtot%pw_w,2)
     235            0 :             do st=1,stars%ng3
     236            0 :                vTot%pw_w(st,js) = vTot%pw_w(st,js) * stars%nstr(st)
     237              :             enddo
     238              :          enddo
     239              :       else
     240            0 :          call juDFT_error("undo vgen_finalize not implemented for noco")
     241              :       endif
     242            0 :    end subroutine undo_vgen_finalize
     243              : 
     244          674 :    subroutine set_kinED(fmpi,   sphhar, atoms, sym,  xcpot, &
     245              :                         input, noco,   stars, vacuum ,cell,     den,     EnergyDen, vTot,kinED)
     246              :       use m_types
     247              :       use m_cdn_io
     248              :       implicit none
     249              :       TYPE(t_mpi),INTENT(IN)       :: fmpi
     250              :       TYPE(t_sphhar),INTENT(IN)    :: sphhar
     251              :       TYPE(t_atoms),INTENT(IN)     :: atoms
     252              :       TYPE(t_sym), INTENT(IN)      :: sym
     253              :       CLASS(t_xcpot),INTENT(IN)    :: xcpot
     254              :       TYPE(t_input),INTENT(IN)     :: input
     255              :       TYPE(t_noco),INTENT(IN)      :: noco
     256              :       TYPE(t_stars),INTENT(IN)     :: stars
     257              :       TYPE(t_vacuum),INTENT(IN)    :: vacuum
     258              :        
     259              :       TYPE(t_cell),INTENT(IN)      :: cell
     260              :       TYPE(t_potden),INTENT(IN)    :: den, EnergyDen, vTot
     261              :       TYPE(t_kinED),INTENT(OUT)    :: kinED
     262              : 
     263          674 :       TYPE(t_potden)               :: vTot_corrected
     264              : 
     265              :       LOGICAL :: perform_MetaGGA
     266          674 :       TYPE(t_potden)    :: core_den, val_den
     267              :       real :: rdum, tempDistance
     268              :       logical :: ldum
     269              :       perform_MetaGGA = ALLOCATED(EnergyDen%mt) &
     270          674 :                       .AND. (xcpot%exc_is_MetaGGA() .or. xcpot%vx_is_MetaGGA())
     271          674 :       if(.not.perform_MetaGGA) return
     272              : #ifdef CPP_LIBXC
     273            0 :       core_den=den
     274            0 :       CALL core_den%resetPotDen
     275              :       call readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym ,&
     276              :                           CDN_ARCHIVE_TYPE_CDN_const,CDN_INPUT_DEN_const,&
     277            0 :                            0,rdum,tempDistance,ldum,core_den,inFilename='cdnc')
     278            0 :       call val_den%subPotDen(den,core_den)
     279              : 
     280            0 :       call vTot_corrected%copyPotDen(vTot)
     281            0 :       call undo_vgen_finalize(vTot_corrected, atoms, noco, stars)
     282              : 
     283            0 :       call set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot_corrected,kinED)
     284              :       call set_kinED_mt(fmpi,   sphhar,    atoms, sym, noco,core_den, val_den, &
     285            0 :                            xcpot, EnergyDen, input, vTot_corrected,kinED)
     286              : #endif
     287          674 :    end subroutine set_kinED
     288              : #ifdef CPP_LIBXC
     289            0 :    subroutine set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot,kinED)
     290              :       use m_types
     291              :       use m_pw_tofrom_grid
     292              :       implicit none
     293              :       CLASS(t_xcpot),INTENT(IN)    :: xcpot
     294              :       TYPE(t_input),INTENT(IN)     :: input
     295              :       TYPE(t_noco),INTENT(IN)      :: noco
     296              :       TYPE(t_stars),INTENT(IN)     :: stars
     297              :       TYPE(t_sym), INTENT(IN)      :: sym
     298              :       TYPE(t_cell),INTENT(IN)      :: cell
     299              :       TYPE(t_potden),INTENT(IN)    :: den, EnergyDen, vTot
     300              :       TYPE(t_kinED),INTENT(INOUT)  :: kinED
     301              : 
     302              :       !local arrays
     303            0 :       REAL, ALLOCATABLE            :: den_rs(:,:), ED_rs(:,:), vTot_rs(:,:)
     304            0 :       TYPE(t_gradients)            :: tmp_grad
     305              : 
     306            0 :       CALL init_pw_grid(stars,sym,cell,xcpot)
     307              : 
     308              :       CALL pw_to_grid(xcpot%needs_grad(), input%jspins, noco%l_noco, stars, &
     309            0 :                       cell,  EnergyDen%pw, tmp_grad, xcpot,    ED_rs)
     310              :       CALL pw_to_grid(xcpot%needs_grad(), input%jspins, noco%l_noco, stars, &
     311            0 :                       cell,  vTot%pw,      tmp_grad, xcpot,    vTot_rs)
     312              :       CALL pw_to_grid(xcpot%needs_grad(), input%jspins, noco%l_noco, stars, &
     313            0 :                       cell,  den%pw,       tmp_grad, xcpot,   den_rs)
     314              : 
     315            0 :       CALL finish_pw_grid()
     316              : 
     317            0 :       call calc_kinEnergyDen_pw(ED_rs, vTot_rs, den_rs, kinED%is)
     318              :       !xcpot%kinED%is  = ED_RS - vTot_RS * den_RS
     319            0 :       kinED%set = .True.
     320            0 :    end subroutine set_kinED_is
     321              : 
     322            0 :    subroutine set_kinED_mt(fmpi,   sphhar,    atoms, sym, noco,core_den, val_den, &
     323              :                            xcpot, EnergyDen, input, vTot,kinED)
     324              :       use m_types
     325              :       use m_mt_tofrom_grid
     326              :       implicit none
     327              :       TYPE(t_mpi),INTENT(IN)         :: fmpi
     328              :       TYPE(t_sphhar),INTENT(IN)      :: sphhar
     329              :       TYPE(t_atoms),INTENT(IN)       :: atoms
     330              :       TYPE(t_sym), INTENT(IN)        :: sym
     331              :       TYPE(t_noco), INTENT(IN)       :: noco
     332              :       TYPE(t_potden),INTENT(IN)      :: core_den, val_den, EnergyDen, vTot
     333              :       CLASS(t_xcpot),INTENT(IN)      :: xcpot
     334              :       TYPE(t_input),INTENT(IN)       :: input
     335              :       TYPE(t_kinED),INTENT(INOUT)    :: kinED
     336              :       INTEGER                        :: jr, loc_n, n, n_start, n_stride, cnt
     337            0 :       REAL,ALLOCATABLE               :: vTot_mt(:,:,:), ED_rs(:,:), vTot_rs(:,:), vTot0_rs(:,:),&
     338              :                                         core_den_rs(:,:), val_den_rs(:,:)
     339            0 :       TYPE(t_gradients)              :: tmp_grad
     340            0 :       TYPE(t_sphhar)                 :: tmp_sphhar
     341              : 
     342              : #ifdef CPP_MPI
     343            0 :       n_start=fmpi%irank+1
     344            0 :       n_stride=fmpi%isize
     345              : #else
     346              :       n_start=1
     347              :       n_stride=1
     348              : #endif
     349            0 :       CALL init_mt_grid(input%jspins,atoms,sphhar,xcpot%needs_grad(),sym)
     350            0 :       loc_n = 0
     351            0 :       allocate(ED_rs(atoms%nsp()*atoms%jmtd, input%jspins))
     352            0 :       allocate(vTot_rs, mold=ED_rs)
     353            0 :       allocate(vTot0_rs, mold=ED_rs)
     354            0 :       allocate(core_den_rs, mold=ED_rs)
     355            0 :       allocate(val_den_rs, mold=ED_rs)
     356              : 
     357              :       call kinED%alloc_mt(atoms%nsp()*atoms%jmtd, input%jspins, &
     358            0 :                                 n_start,                atoms%ntype,  n_stride)
     359            0 :       loc_n = 0
     360            0 :       do n = n_start,atoms%ntype,n_stride
     361            0 :          loc_n = loc_n + 1
     362              : 
     363            0 :          if(.not. allocated(vTot_mt)) then
     364              :             allocate(vTot_mt(lbound(vTot%mt, dim=1):ubound(vTot%mt, dim=1),&
     365              :                              lbound(vTot%mt, dim=2):ubound(vTot%mt, dim=2),&
     366            0 :                              lbound(vTot%mt, dim=4):ubound(vTot%mt, dim=4)))
     367              :          endif
     368              : 
     369            0 :          do jr=1,atoms%jri(n)
     370            0 :             vTot_mt(jr,0:,:) = vTot%mt(jr,0:,n,:) * atoms%rmsh(jr,n)**2
     371              :          enddo
     372              :          CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms,sym, sphhar,.TRUE., EnergyDen%mt(:, 0:, n, :), &
     373            0 :                          n,  noco,   tmp_grad,     ED_rs)
     374              :          CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,sphhar,.TRUE., vTot_mt(:,0:,:), &
     375            0 :                          n,     noco,tmp_grad,     vTot_rs)
     376              : 
     377            0 :          tmp_sphhar%nlhd = sphhar%nlhd
     378            0 :          tmp_sphhar%nlh  = [(0, cnt=1,size(sphhar%nlh))]
     379              : 
     380              :          CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,tmp_sphhar,.TRUE., vTot_mt(:,0:0,:), &
     381            0 :                          n,    noco, tmp_grad,     vTot0_rs)
     382              :          CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,sphhar,.TRUE., &
     383            0 :                          core_den%mt(:,0:,n,:), n,noco, tmp_grad, core_den_rs)
     384              :          CALL mt_to_grid(xcpot%needs_grad(), input%jspins, atoms, sym,sphhar,.TRUE., &
     385            0 :                          val_den%mt(:,0:,n,:), n,noco, tmp_grad, val_den_rs)
     386              : 
     387              :          call calc_kinEnergyDen_mt(ED_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, &
     388            0 :                                    kinED%mt(:,:,loc_n))
     389              :          !xcpot%kinED%mt(:,:,loc_n) = ED_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
     390              :       enddo
     391            0 :       kinED%set = .True.
     392            0 :       CALL finish_mt_grid()
     393            0 :    end subroutine set_kinED_mt
     394              : #endif
     395            0 : END MODULE m_metagga
        

Generated by: LCOV version 2.0-1