LCOV - code coverage report
Current view: top level - xc-pot - metagga.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 100 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 10 0.0 %

          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             : CONTAINS
      17           0 :    SUBROUTINE calc_kinEnergyDen_pw(EnergyDen_rs, vTot_rs, den_rs, kinEnergyDen_RS)
      18             :       USE m_juDFT_stop
      19             :       !use m_cdngen
      20             :       IMPLICIT NONE
      21             :       REAL, INTENT(in)                 :: den_RS(:,:), EnergyDen_RS(:,:), vTot_RS(:,:)
      22             :       REAL, INTENT(inout), allocatable :: kinEnergyDen_RS(:,:)
      23             : #ifdef CPP_LIBXC
      24             :       REAL, PARAMETER                  :: eps = 1e-15
      25             : 
      26           0 :       kinEnergyDen_RS = EnergyDen_RS - vTot_RS * den_RS
      27             : #else
      28             :       CALL juDFT_error("MetaGGA require LibXC",hint="compile Fleur with LibXC (e.g. by giving '-external libxc' to ./configure")
      29             : #endif
      30           0 :    END SUBROUTINE calc_kinEnergyDen_pw
      31             : 
      32           0 :    SUBROUTINE calc_kinEnergyDen_mt(EnergyDen_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, &
      33           0 :                                    kinEnergyDen_RS)
      34             :       USE m_juDFT_stop
      35             :       USE m_juDFT_string
      36             :       implicit none
      37             :       REAL, INTENT(in)                 :: EnergyDen_RS(:,:), vTot_rs(:,:), vTot0_rs(:,:), core_den_rs(:,:), val_den_rs(:,:)
      38             :       REAL, INTENT(inout)              :: kinEnergyDen_RS(:,:)
      39             : 
      40             : #ifdef CPP_LIBXC
      41           0 :       kinEnergyDen_RS = EnergyDen_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
      42             : #else
      43             :       CALL juDFT_error("MetaGGA require LibXC",hint="compile Fleur with LibXC (e.g. by giving '-external libxc' to ./configure")
      44             : #endif
      45           0 :    END SUBROUTINE calc_kinEnergyDen_mt
      46             : 
      47             : 
      48           0 :    SUBROUTINE calc_EnergyDen(eig_id, mpi, kpts, noco, input, banddos, cell, atoms, enpara, stars, &
      49             :          vacuum, DIMENSION, sphhar, sym, vTot, oneD, results, EnergyDen)
      50             :       ! calculates the energy density
      51             :       ! EnergyDen = \sum_i n_i(r) \varepsilon_i
      52             :       ! where n_i(r) is the one-particle density
      53             :       ! and \varepsilon_i are the eigenenergies
      54             : 
      55             :       USE m_types_setup
      56             :       USE m_types_potden
      57             :       USE m_types_kpts
      58             :       USE m_types_mpi
      59             :       USE m_types_enpara
      60             :       USE m_types_misc
      61             :       USE m_types_regionCharges
      62             :       USE m_types_dos
      63             :       USE m_types_cdnval
      64             :       USE m_cdnval
      65             : 
      66             :       IMPLICIT NONE
      67             : 
      68             :       INTEGER,           INTENT(in)           :: eig_id
      69             :       TYPE(t_mpi),       INTENT(in)           :: mpi
      70             :       TYPE(t_kpts),      INTENT(in)           :: kpts
      71             :       TYPE(t_noco),      INTENT(in)           :: noco
      72             :       TYPE(t_input),     INTENT(in)           :: input
      73             :       TYPE(t_banddos),   INTENT(in)           :: banddos
      74             :       TYPE(t_cell),      INTENT(in)           :: cell
      75             :       TYPE(t_atoms),     INTENT(in)           :: atoms
      76             :       TYPE(t_enpara),    INTENT(in)           :: enpara
      77             :       TYPE(t_stars),     INTENT(in)           :: stars
      78             :       TYPE(t_vacuum),    INTENT(in)           :: vacuum
      79             :       TYPE(t_dimension), INTENT(in)           :: DIMENSION
      80             :       TYPE(t_sphhar),    INTENT(in)           :: sphhar
      81             :       TYPE(t_sym),       INTENT(in)           :: sym
      82             :       TYPE(t_potden),    INTENT(in)           :: vTot
      83             :       TYPE(t_oneD),      INTENT(in)           :: oneD
      84             :       TYPE(t_results),   INTENT(in)           :: results
      85             :       TYPE(t_potden),    INTENT(inout)        :: EnergyDen
      86             : 
      87             :       ! local
      88             :       INTEGER                         :: jspin
      89             : 
      90           0 :       TYPE(t_regionCharges)           :: regCharges
      91           0 :       TYPE(t_dos)                     :: dos
      92           0 :       TYPE(t_moments)                 :: moments
      93           0 :       TYPE(t_results)                 :: tmp_results
      94           0 :       TYPE(t_cdnvalJob)               :: cdnvalJob
      95             :       TYPE(t_potden)                  :: aux_den, real_den
      96             : 
      97             : 
      98           0 :       CALL regCharges%init(input, atoms)
      99           0 :       CALL dos%init(input,        atoms, DIMENSION, kpts, vacuum)
     100           0 :       CALL moments%init(input,    atoms)
     101           0 :       tmp_results = results
     102             : 
     103           0 :       DO jspin = 1,input%jspins
     104           0 :          CALL cdnvalJob%init(mpi,input,kpts,noco,results,jspin)
     105             : 
     106             : 
     107             :          ! replace brillouin weights with auxillary weights
     108           0 :          CALL calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, cdnvalJob%weights)
     109             : 
     110             :          CALL cdnval(eig_id, mpi, kpts, jspin, noco, input, banddos, cell, atoms, &
     111             :             enpara, stars, vacuum, DIMENSION, sphhar, sym, vTot, oneD, cdnvalJob, &
     112           0 :             EnergyDen, regCharges, dos, tmp_results, moments)
     113             :       ENDDO
     114             : 
     115           0 :    END SUBROUTINE calc_EnergyDen
     116             : 
     117           0 :    SUBROUTINE calc_EnergyDen_auxillary_weights(eig_id, kpts, jspin, f_ik)
     118             :       USE m_types_kpts
     119             :       USE m_eig66_io
     120             :       IMPLICIT NONE
     121             :       ! calculates new (auxillary-)weights as
     122             :       ! f_iks = w_iks * E_iks
     123             :       !, where  f_iks are the new (auxillary-)weights
     124             :       ! w_iks are the weights used in brillouin zone integration
     125             :       ! E_iks are the eigen energies
     126             : 
     127             :       INTEGER,      INTENT(in)        :: eig_id
     128             :       INTEGER,      INTENT(in)        :: jspin
     129             :       TYPE(t_kpts), INTENT(in)        :: kpts
     130             :       REAL,         INTENT(inout)     :: f_ik(:,:) ! f_ik(band_idx, kpt_idx)
     131             : 
     132             :       ! local vars
     133           0 :       REAL                       :: e_i(SIZE(f_ik,dim=1))
     134             :       INTEGER                    :: ikpt
     135             : 
     136           0 :       DO ikpt = 1,kpts%nkpt
     137           0 :          CALL read_eig(eig_id,ikpt,jspin, eig=e_i)
     138           0 :          f_ik(:,ikpt) = f_ik(:,ikpt) * e_i
     139             :       ENDDO
     140           0 :    END SUBROUTINE calc_EnergyDen_auxillary_weights
     141             : 
     142           0 :    subroutine set_zPrime(dim_idx, zMat, kpt, lapw, cell, zPrime)
     143             :       USE m_types
     144             :       USE m_constants
     145             :       implicit none
     146             :       INTEGER, intent(in)      :: dim_idx
     147             :       TYPE (t_mat), intent(in) :: zMat
     148             :       REAL, intent(in)         :: kpt(3) 
     149             :       TYPE(t_lapw), intent(in) :: lapw
     150             :       TYPE(t_cell), intent(in) :: cell
     151             :       TYPE (t_mat)             :: zPrime
     152             : 
     153             :       REAL                     :: k_plus_g(3), fac
     154             :       INTEGER                  :: basis_idx
     155             : 
     156           0 :       call zPrime%free()
     157           0 :       call zPrime%init(zMat)
     158             : 
     159           0 :       do basis_idx = 1,size(lapw%gvec,dim=2)
     160           0 :          k_plus_g = kpt + lapw%gvec(:,basis_idx,1)
     161           0 :          k_plus_g = internal_to_rez(cell, k_plus_g)
     162             : 
     163           0 :          fac = k_plus_g(dim_idx)
     164           0 :          if(zPrime%l_real) then
     165           0 :             zPrime%data_r(basis_idx,:) =            fac * zMat%data_r(basis_idx,:) 
     166             :          else
     167           0 :             zPrime%data_c(basis_idx,:) = ImagUnit * fac * zMat%data_c(basis_idx,:) 
     168             :          endif
     169             :       enddo
     170           0 :    end subroutine set_zPrime
     171             : 
     172           0 :    function internal_to_rez(cell, vec) result(res)
     173             :       use m_types
     174             :       implicit none
     175             :       type(t_cell), intent(in) :: cell
     176             :       real, intent(in)      :: vec(3)
     177             :       real                  :: res(3)
     178             : 
     179           0 :       res = matmul(transpose(cell%bmat), vec)
     180           0 :    end function internal_to_rez
     181             : 
     182           0 :    subroutine undo_vgen_finalize(vtot, atoms, noco, stars)
     183             :       use m_types
     184             :       use m_constants
     185             :       use m_judft
     186             :       implicit none
     187             :       TYPE(t_potden), intent(inout)  :: vtot
     188             :       type(t_atoms), intent(in)      :: atoms
     189             :       type(t_noco), intent(in)       :: noco
     190             :       type(t_stars), intent(in)      :: stars
     191             :    
     192             :       integer                        :: js, n, st
     193             : 
     194           0 :       do js = 1,size(vtot%mt,4)
     195           0 :          do n = 1,atoms%ntype
     196             :             vTot%mt(:atoms%jri(n),0,n,js) = vtot%mt(:atoms%jri(n),0,n,js) &
     197           0 :                   / (atoms%rmsh(:atoms%jri(n),n) / sfp_const )
     198             :          enddo
     199             :       enddo
     200             : 
     201           0 :       if(.not. noco%l_noco) then
     202           0 :          do js=1,size(vtot%pw_w,2)
     203           0 :             do st=1,stars%ng3
     204           0 :                vTot%pw_w(st,js) = vTot%pw_w(st,js) * stars%nstr(st)
     205             :             enddo
     206             :          enddo
     207             :       else
     208           0 :          call juDFT_error("undo vgen_finalize not implemented for noco")
     209             :       endif
     210           0 :    end subroutine undo_vgen_finalize
     211             : 
     212           0 :    subroutine set_kinED(mpi,   sphhar, atoms, sym, core_den, val_den, xcpot, &
     213             :                         input, noco,   stars, cell,     den,     EnergyDen, vTot)
     214             :       use m_types
     215             :       implicit none
     216             :       TYPE(t_mpi),INTENT(IN)       :: mpi
     217             :       TYPE(t_sphhar),INTENT(IN)    :: sphhar
     218             :       TYPE(t_atoms),INTENT(IN)     :: atoms
     219             :       TYPE(t_sym), INTENT(IN)      :: sym
     220             :       TYPE(t_potden),INTENT(IN)    :: core_den, val_den
     221             :       CLASS(t_xcpot),INTENT(INOUT) :: xcpot
     222             :       TYPE(t_input),INTENT(IN)     :: input
     223             :       TYPE(t_noco),INTENT(IN)      :: noco
     224             :       TYPE(t_stars),INTENT(IN)     :: stars
     225             :       TYPE(t_cell),INTENT(IN)      :: cell
     226             :       TYPE(t_potden),INTENT(IN)    :: den, EnergyDen, vTot
     227             :       
     228           0 :       TYPE(t_potden)               :: vTot_corrected
     229             :    
     230           0 :       call vTot_corrected%copyPotDen(vTot)
     231           0 :       call undo_vgen_finalize(vTot_corrected, atoms, noco, stars)
     232             : 
     233           0 :       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           0 :                            xcpot, EnergyDen, input, vTot_corrected)
     236           0 :    end subroutine set_kinED
     237             : 
     238           0 :    subroutine set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot)
     239             :       use m_types
     240             :       use m_pw_tofrom_grid
     241             :       implicit none
     242             :       CLASS(t_xcpot),INTENT(INOUT) :: xcpot
     243             :       TYPE(t_input),INTENT(IN)     :: input
     244             :       TYPE(t_noco),INTENT(IN)      :: noco
     245             :       TYPE(t_stars),INTENT(IN)     :: stars
     246             :       TYPE(t_sym), INTENT(IN)      :: sym
     247             :       TYPE(t_cell),INTENT(IN)      :: cell
     248             :       TYPE(t_potden),INTENT(IN)    :: den, EnergyDen, vTot
     249             : 
     250             :       !local arrays
     251           0 :       REAL, ALLOCATABLE            :: den_rs(:,:), ED_rs(:,:), vTot_rs(:,:)
     252           0 :       TYPE(t_gradients)            :: tmp_grad
     253             :       
     254           0 :       CALL init_pw_grid(xcpot,stars,sym,cell)
     255             : 
     256             :       CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
     257             :                       cell,  EnergyDen%pw, tmp_grad,    ED_rs)
     258             :       CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
     259             :                       cell,  vTot%pw,      tmp_grad,    vTot_rs)
     260             :       CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
     261             :                       cell,  den%pw,       tmp_grad,    den_rs)
     262             : 
     263           0 :       CALL finish_pw_grid()
     264             :       
     265           0 :       call calc_kinEnergyDen_pw(ED_rs, vTot_rs, den_rs, xcpot%kinED%is)
     266             :       !xcpot%kinED%is  = ED_RS - vTot_RS * den_RS
     267           0 :       xcpot%kinED%set = .True.
     268           0 :    end subroutine set_kinED_is
     269             : 
     270           0 :    subroutine set_kinED_mt(mpi,   sphhar,    atoms, sym, core_den, val_den, &
     271             :                            xcpot, EnergyDen, input, vTot)
     272             :       use m_types
     273             :       use m_mt_tofrom_grid
     274             :       implicit none
     275             :       TYPE(t_mpi),INTENT(IN)         :: mpi
     276             :       TYPE(t_sphhar),INTENT(IN)      :: sphhar
     277             :       TYPE(t_atoms),INTENT(IN)       :: atoms
     278             :       TYPE(t_sym), INTENT(IN)        :: sym
     279             :       TYPE(t_potden),INTENT(IN)      :: core_den, val_den, EnergyDen, vTot
     280             :       CLASS(t_xcpot),INTENT(INOUT)   :: xcpot
     281             :       TYPE(t_input),INTENT(IN)       :: input
     282             : 
     283             :       INTEGER                        :: jr, loc_n, n, n_start, n_stride, cnt
     284           0 :       REAL,ALLOCATABLE               :: vTot_mt(:,:,:), ED_rs(:,:), vTot_rs(:,:), vTot0_rs(:,:),&
     285           0 :                                         core_den_rs(:,:), val_den_rs(:,:)
     286           0 :       TYPE(t_gradients)              :: tmp_grad
     287           0 :       TYPE(t_sphhar)                 :: tmp_sphhar
     288             : 
     289             : #ifdef CPP_MPI
     290           0 :       n_start=mpi%irank+1
     291           0 :       n_stride=mpi%isize
     292             : #else
     293             :       n_start=1
     294             :       n_stride=1
     295             : #endif
     296           0 :       CALL init_mt_grid(input%jspins,atoms,sphhar,xcpot,sym)
     297           0 :       loc_n = 0
     298           0 :       allocate(ED_rs(atoms%nsp()*atoms%jmtd, input%jspins))
     299           0 :       allocate(vTot_rs, mold=ED_rs)
     300           0 :       allocate(vTot0_rs, mold=ED_rs)
     301           0 :       allocate(core_den_rs, mold=ED_rs)
     302           0 :       allocate(val_den_rs, mold=ED_rs)
     303             : 
     304             :       call xcpot%kinED%alloc_mt(atoms%nsp()*atoms%jmtd, input%jspins, &
     305           0 :                                 n_start,                atoms%ntype,  n_stride)
     306           0 :       loc_n = 0
     307           0 :       do n = n_start,atoms%ntype,n_stride
     308           0 :          loc_n = loc_n + 1
     309             : 
     310           0 :          if(.not. allocated(vTot_mt)) then
     311             :             allocate(vTot_mt(lbound(vTot%mt, dim=1):ubound(vTot%mt, dim=1),&
     312             :                              lbound(vTot%mt, dim=2):ubound(vTot%mt, dim=2),&
     313           0 :                              lbound(vTot%mt, dim=4):ubound(vTot%mt, dim=4)))
     314             :          endif
     315             :          
     316           0 :          do jr=1,atoms%jri(n)
     317           0 :             vTot_mt(jr,0:,:) = vTot%mt(jr,0:,n,:) * atoms%rmsh(jr,n)**2
     318             :          enddo
     319             :          CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, EnergyDen%mt(:, 0:, n, :), &
     320           0 :                          n,     tmp_grad,     ED_rs)
     321             :          CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, vTot_mt(:,0:,:), &
     322           0 :                          n,     tmp_grad,     vTot_rs)
     323             :          
     324           0 :          tmp_sphhar%nlhd = sphhar%nlhd
     325           0 :          tmp_sphhar%nlh  = [(0, cnt=1,size(sphhar%nlh))]
     326             : 
     327             :          CALL mt_to_grid(xcpot, input%jspins, atoms, tmp_sphhar, vTot_mt(:,0:0,:), &
     328           0 :                          n,     tmp_grad,     vTot0_rs)
     329             :          CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, &
     330           0 :                          core_den%mt(:,0:,n,:), n, tmp_grad, core_den_rs)
     331             :          CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, &
     332           0 :                          val_den%mt(:,0:,n,:), n, tmp_grad, val_den_rs)
     333             :          
     334             :          call calc_kinEnergyDen_mt(ED_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, &
     335           0 :                                    xcpot%kinED%mt(:,:,loc_n))
     336             :          !xcpot%kinED%mt(:,:,loc_n) = ED_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs)
     337             :       enddo
     338           0 :       xcpot%kinED%set = .True.
     339           0 :       CALL finish_mt_grid()
     340           0 :    end subroutine set_kinED_mt
     341             : END MODULE m_metagga

Generated by: LCOV version 1.13