LCOV - code coverage report
Current view: top level - src/67_common - m_xc_tb09.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 86.4 % 169 146
Test Date: 2026-09-19 15:24:51 Functions: 100.0 % 1 1

            Line data    Source code
       1              : !!****m* m_xc_tb09/m_xc_tb09
       2              : !! NAME
       3              : !!  m_xc_tb09
       4              : !!
       5              : !! FUNCTION
       6              : !!  This module contains routine(s) related to exchange-correlation Tran-Blaha 2009
       7              : !!    functional (modified Becke-Johnson)
       8              : !!
       9              : !! COPYRIGHT
      10              : !! Copyright (C) 2023-2026 ABINIT group (MT)
      11              : !! This file is distributed under the terms of the
      12              : !! GNU General Public License, see ~abinit/COPYING
      13              : !! or http://www.gnu.org/copyleft/gpl.txt .
      14              : !!
      15              : !! SOURCE
      16              : 
      17              : #if defined HAVE_CONFIG_H
      18              : #include "config.h"
      19              : #endif
      20              : 
      21              : #include "abi_common.h"
      22              : 
      23              : MODULE m_xc_tb09
      24              : 
      25              :  use defs_basis
      26              :  use m_abicore
      27              :  use m_xmpi
      28              :  use m_errors
      29              :  use m_time,         only : timab
      30              :  use defs_abitypes,  only : MPI_type
      31              : 
      32              :  use m_xctk,         only : xcden
      33              :  use m_drivexc,      only : mkdenpos
      34              :  use m_geometry,     only : metric
      35              : 
      36              :  use m_pawang,       only : pawang_type
      37              :  use m_pawrad,       only : pawrad_type,pawrad_deducer0,nderiv_gen,simp_gen
      38              :  use m_pawtab,       only : pawtab_type
      39              :  use m_pawrhoij,     only : pawrhoij_type
      40              :  use m_paw_denpot,   only : pawdensities
      41              :  use m_paral_atom,   only : get_my_atmtab,free_my_atmtab
      42              : 
      43              :  use libxc_functionals, only : libxc_functional_type,libxc_functionals_set_c_tb09, &
      44              : &                              libxc_functionals_is_tb09,libxc_functionals_ixc
      45              : 
      46              :  implicit none
      47              : 
      48              :  private
      49              : 
      50              : !public procedure(s)
      51              :  public :: xc_tb09_update_c ! Compute the C parameter for the TB09 functional (NCPP+PAW)
      52              : 
      53              : !Private variables
      54              :  real(dp),save :: grho_over_rho_pw    ! Contribution to TB09 c from the plane-wave pseudo-density
      55              :  real(dp),save :: grho_over_rho_paw   ! Contribution to TB09 c from the on-site PAW densities
      56              : 
      57              : CONTAINS  !========================================================================================
      58              : !!***
      59              : 
      60              : !----------------------------------------------------------------------
      61              : 
      62              : !!****f* m_xc_tb09/xc_tb09_update_c
      63              : !! NAME
      64              : !! xc_tb09_update_c
      65              : !!
      66              : !! FUNCTION
      67              : !! This routine computes from the density the C parameter of the Tran-Blaha 2009
      68              : !!  XC functional (modified Becke-Johnson)
      69              : !!  C = 1/V Int[|Grad(Rho)|/Rho] + Sum_PAW_aug_regions[ Natom/V [|Grad(Rho)|/Rho] ]
      70              : !!
      71              : !! INPUTS
      72              : !!  intxc = 1 if the XC functional has to be interpolated on a more refined mesh
      73              : !!  ixc= choice of exchange-correlation scheme
      74              : !!  mpi_enreg= information about MPI parallelization
      75              : !!  natom= total number of atoms in cell
      76              : !!  nfft= (effective) number of FFT grid points (for this processor)
      77              : !!  ngfft(18)= contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
      78              : !!  nhat(nfft,nspden*nhatdim)= -PAW only- compensation density
      79              : !!  nhatdim= -PAW only- 0 if nhat array is not used ; 1 otherwise
      80              : !!  nhatgr(nfft,xcdata%nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
      81              : !!  nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
      82              : !!  nspden= number of spin-density components
      83              : !!  ntypat= number of types of atoms in unit cell.
      84              : !!  n3xccc= dimension of the xccc3d array (0 or nfft or cplx*nfft).
      85              : !!  pawang <type(pawang_type)> =paw angular mesh and related data
      86              : !!  pawrad(ntypat) <type(pawrad_type)> =paw radial mesh and related data
      87              : !!  pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
      88              : !!  pawtab(ntypat) <type(pawtab_type)>= paw tabulated starting data
      89              : !!  pawxcdev= choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
      90              : !!  rhor(nfft,nspden)= electron density in real space in electrons/bohr**3
      91              : !!  rprimd(3,3)=dimensional primitive translations in real space (bohr)
      92              : !!  usepaw = 1 for PAW computation, 0 else
      93              : !!  xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
      94              : !!  xc_denpos= lowest allowed density
      95              : !!
      96              : !!  ----- Optional arguments -----
      97              : !!  [computation_type]= 'only_pw' : update only plane-wave contribution
      98              : !!                      'only_paw': update only PAW on-site contribution
      99              : !!                      'all' (default): compute all contributions
     100              : !!  [mpi_atmtab(:)]= indexes of the atoms treated by current proc
     101              : !!  [comm_atom]= MPI communicator over atoms
     102              : !!  [xc_funcs(2)]= <type(libxc_functional_type)>=libxc XC functionals (if not the global one)
     103              : !!
     104              : !! OUTPUT
     105              : !!  No output
     106              : !!  The C parameter is directly set in the xc_functional datastructure
     107              : !!
     108              : !! SOURCE
     109              : 
     110       110324 : subroutine xc_tb09_update_c(intxc,ixc,mpi_enreg,natom,nfft,ngfft,nhat,nhatdim, &
     111       110324 : &                           nhatgr,nhatgrdim,nspden,ntypat,n3xccc,pawang,pawrad, &
     112        55162 : &                           pawrhoij,pawtab,pawxcdev,rhor,rprimd,usepaw,xccc3d,xc_denpos, &
     113        55162 : &                           computation_type,mpi_atmtab,comm_atom,xc_funcs) ! Optional arguments
     114              : 
     115              : !Arguments ------------------------------------
     116              : !scalars
     117              :  integer,intent(in) :: intxc,ixc,n3xccc,nfft,natom,nhatdim,nhatgrdim
     118              :  integer,intent(in) :: nspden,ntypat,pawxcdev,usepaw
     119              :  integer,intent(in),optional :: comm_atom
     120              :  real(dp),intent(in),optional :: xc_denpos
     121              :  character(len=*),intent(in),optional :: computation_type
     122              :  type(pawang_type),intent(in) :: pawang
     123              :  type(MPI_type),intent(in) :: mpi_enreg
     124              : !arrays
     125              :  integer,intent(in) :: ngfft(18)
     126              :  integer,intent(in),optional,target :: mpi_atmtab(:)
     127              :  real(dp),intent(in) :: nhat(nfft,nspden*nhatdim),nhatgr(nfft,nspden,3*nhatgrdim)
     128              :  real(dp),intent(in) :: rhor(nfft,nspden)
     129              :  real(dp),intent(in) :: rprimd(3,3),xccc3d(n3xccc)
     130              :  type(libxc_functional_type),intent(inout),optional :: xc_funcs(2)
     131              :  type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw)
     132              :  type(pawrhoij_type),intent(in) :: pawrhoij(:)
     133              :  type(pawtab_type),intent(in),target :: pawtab(ntypat*usepaw)
     134              : 
     135              : !Local variables-------------------------------
     136              : !scalars
     137              :  integer :: cplex,iatom,iatom_tot,ierr,ii,ilm,iloop,ir,itypat,ixc_from_lib,ipts
     138              :  integer :: ishift,iwarn,lm_size,lm_size_eff,mesh_size,my_comm_atom,my_natom,ngrad
     139              :  integer :: nfftot,npts,nspden1,nzlmopt,option,opt_compch,opt_dens,pawprtvol,usecore
     140              :  integer :: usenhat,usexcnhat
     141              :  real(dp) :: compch_dum,factor,grho_over_rho,grho2,rho,sumg,ucvol,xc_c_tb09
     142              :  logical :: my_atmtab_allocated,paral_atom,test_tb09,use_exact_nhat_gradient
     143              :  logical :: compute_pw,compute_paw
     144              :  character(len=500) :: msg
     145              : !arrays
     146        55162 :  integer,pointer :: my_atmtab(:)
     147        55162 :  logical,allocatable :: lmselect(:),lmselect_tmp(:)
     148              :  real(dp) :: gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3)
     149              : !real(dp) :: tsec(2)
     150        55162 :  real(dp),allocatable :: dylmdr(:,:,:),ff(:)
     151        55162 :  real(dp),allocatable :: rhotot(:),grhotot(:,:),rhonow(:,:,:)
     152        55162 :  real(dp),allocatable :: rhoxc(:,:),drhoxc(:),dcorexc(:)
     153        55162 :  real(dp),allocatable :: rho1(:,:,:),trho1(:,:,:),nhat1(:,:,:)
     154        55162 :  real(dp), ABI_CONTIGUOUS pointer :: corexc(:)
     155              : 
     156              : ! *************************************************************************
     157              : 
     158              : !call timab(xx,1,tsec)
     159              : 
     160              : !Type of calculation
     161        55162 :  compute_pw=.true. ; compute_paw=.true.
     162        55162 :  if (present(computation_type)) then
     163        54446 :    if (computation_type=='only_pw') compute_paw=.false.
     164        54446 :    if (computation_type=='only_paw') compute_pw=.false.
     165              :  end if
     166            0 :  if ((.not.compute_pw).and.(.not.compute_paw)) return
     167              : 
     168              : !Nothing to do if XC functional is not TB09
     169        55162 :  if (present(xc_funcs)) then
     170            0 :    test_tb09=libxc_functionals_is_tb09(xc_functionals=xc_funcs)
     171              :  else
     172        55162 :    test_tb09=libxc_functionals_is_tb09()
     173              :  end if
     174        55162 :  if (.not.test_tb09) return
     175              : 
     176              : !Check options
     177           71 :  if(ixc<0) then
     178           71 :    if (present(xc_funcs)) then
     179            0 :      ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs)
     180              :    else
     181           71 :      ixc_from_lib=libxc_functionals_ixc()
     182              :    end if
     183           71 :    if (ixc /= ixc_from_lib) then
     184              :      write(msg, '(a,i0,a,a,i0)')&
     185            0 : &     'The value of ixc specified in input, ixc = ',ixc,ch10,&
     186            0 : &     'differs from the one used to initialize the functional ',ixc_from_lib
     187            0 :      ABI_BUG(msg)
     188              :    end if
     189              :  end if
     190           71 :  if (usepaw==1.and.compute_paw) then
     191           31 :    if (pawxcdev/=0) then
     192            0 :      msg='TB09 (mBJ) XC functional only available with pawxcdev=0!'
     193            0 :      ABI_BUG(msg)
     194              :    end if
     195           31 :    if(pawang%angl_size==0) then
     196            0 :      msg='pawang%angl_size=0!'
     197            0 :      ABI_BUG(msg)
     198              :    end if
     199           31 :    if(.not.allocated(pawang%ylmr)) then
     200            0 :      msg='pawang%ylmr must be allocated!'
     201            0 :      ABI_BUG(msg)
     202              :    end if
     203           31 :    if(.not.allocated(pawang%ylmrgr)) then
     204            0 :      msg='pawang%ylmrgr must be allocated!'
     205            0 :      ABI_BUG(msg)
     206              :    end if
     207           31 :    if (size(pawrhoij)/=mpi_enreg%my_natom) then
     208            0 :      msg='Invalid pawrhoij size!'
     209            0 :      ABI_BUG(msg)
     210              :    end if
     211              :  end if
     212              : 
     213              : !Initialize cell data
     214           71 :  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
     215              : 
     216              : !Other initializations
     217          102 :  usexcnhat=0; if (usepaw==1) usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
     218              : 
     219              : !Initialization of Int[|Grad(Rho)|/Rho] contributions (plane-wave and on-site if PAW)
     220           71 :  if (compute_pw) grho_over_rho_pw=zero
     221           71 :  if (compute_paw) grho_over_rho_paw=zero
     222              : 
     223              : ! ========================================================================
     224              : ! =========== Plane wave contribution to Int[|Grad(Rho)|/Rho] ============
     225              : 
     226           71 :  if (compute_pw) then
     227              : 
     228              : !  Total density used in XC (with/without core density and compensation density)
     229          213 :    ABI_MALLOC(rhotot,(nfft))
     230      1394031 :    rhotot(1:nfft)=rhor(1:nfft,1)
     231       428615 :    if (usexcnhat==0.and.nhatdim==1) rhotot(1:nfft)=rhotot(1:nfft)-nhat(1:nfft,1)
     232       428615 :    if (n3xccc>0) rhotot(1:nfft)=rhotot(1:nfft)+xccc3d(1:nfft)
     233              : 
     234              : !  When possible, it is better to use the exact (analytical) gradient of nhat
     235              : !  In that case, we substract nhat from the density now
     236              : !  If not possible, the gradient will be computed in Fourier space (see call to xcden)
     237           71 :    use_exact_nhat_gradient=(nhatdim==1.and.nhatgrdim==1.and.usexcnhat==1.and.intxc==0)
     238            0 :    if (use_exact_nhat_gradient) rhotot(1:nfft)=rhotot(1:nfft)-nhat(1:nfft,1)
     239              : 
     240              : !  Loop on unshifted or shifted grids
     241          142 :    do ishift=0,intxc
     242              : 
     243              : !    Set up density on unshifted or shifted grid (will be in rhonow(:,:,1)),
     244              : !    as well as the gradient of the density, also on the unshifted
     245              : !    or shifted grid (will be in rhonow(:,:,2:4)), if needed.
     246          213 :      ABI_MALLOC(rhonow,(nfft,1,4))
     247           71 :      nspden1=1 ; cplex=1 ; ngrad=2 ; qphon(:)=zero
     248           71 :      call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden1,qphon,rhotot,rhonow)
     249              : 
     250              : !    PAW: if they were previously substracted, re-add nhat and its "exact" gradient
     251           71 :      if (use_exact_nhat_gradient) then
     252            0 :        rhonow(1:nfft,1,1)=rhonow(1:nfft,1,1)+nhat(1:nfft,1)
     253            0 :        do ii=1,3
     254            0 :          rhonow(1:nfft,1,ii+1)=rhonow(1:nfft,1,ii+1)+nhatgr(1:nfft,1,ii)
     255              :        end do
     256              :      end if
     257              : 
     258              : !    Make the density positive everywhere (but do not care about gradients)
     259           71 :      nspden1=1 ; iwarn=1  ; option=1
     260           71 :      call mkdenpos(iwarn,nfft,nspden1,option,rhonow(:,1,1),xc_denpos)
     261              : 
     262              : !    Accumulate integral of |Grad_rho|/Rho (to be used for TB09 XC)
     263      1394031 :      do ipts=1,nfft
     264      1393960 :        rho=rhonow(ipts,1,1)
     265      1394031 :        if (abs(rho)>tol10) then
     266      1393960 :          grho2=rhonow(ipts,1,2)**2+rhonow(ipts,1,3)**2+rhonow(ipts,1,4)**2
     267      1393960 :          grho_over_rho_pw=grho_over_rho_pw+sqrt(grho2)/rho
     268              :        end if
     269              :      end do
     270              : 
     271          142 :      ABI_FREE(rhonow)
     272              : 
     273              : !  End loop on unshifted or shifted grids
     274              :    end do
     275              : 
     276              : !  Normalize the integral
     277           71 :    nfftot=ngfft(1)*ngfft(2)*ngfft(3)
     278           71 :    grho_over_rho_pw=grho_over_rho_pw*ucvol/dble(nfftot)/dble(intxc+1)
     279              : 
     280              : !  Reduction in case of FFT distribution
     281           71 :    if (mpi_enreg%nproc_fft>1) then
     282            0 :      call xmpi_sum(grho_over_rho_pw,mpi_enreg%comm_fft,ierr)
     283              :    end if
     284              : 
     285           71 :    ABI_FREE(rhotot)
     286              : 
     287              :  end if ! compute_pw
     288              : 
     289              : ! ========================================================================
     290              : ! ========== PAW on-site contributions to Int[|Grad(Rho)|/Rho] ===========
     291              : 
     292           71 :  if (usepaw==1.and.compute_paw) then
     293              : 
     294              : !  Set up parallelism over atoms
     295           31 :    my_natom=mpi_enreg%my_natom
     296           31 :    paral_atom=(present(comm_atom).and.(my_natom/=natom))
     297           31 :    nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
     298           31 :    my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
     299           31 :    call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
     300              : 
     301           31 :    usecore=1 ; nzlmopt=-1
     302           31 :    npts=pawang%angl_size
     303              : 
     304              : !  Compute Ylm gradients
     305              : !    dYlm/dr_i = { dYlm/dr_i^hat - r_i^hat * Sum_j[dYlm/dr_j^hat r_j^hat] } * (1/r)
     306          124 :    ABI_MALLOC(dylmdr,(3,npts,pawang%ylm_size))
     307          527 :    do ilm=1,pawang%ylm_size
     308        77903 :      do ipts=1,npts
     309       309504 :        factor=sum(pawang%ylmrgr(1:3,ilm,ipts)*pawang%anginit(1:3,ipts))
     310       310000 :        dylmdr(1:3,ipts,ilm)=pawang%ylmrgr(1:3,ilm,ipts)-factor*pawang%anginit(1:3,ipts)
     311              :      end do
     312              :    end do
     313              : 
     314              : !  Loop on atom sites
     315           93 :    do iatom=1,my_natom
     316           62 :      iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
     317              : 
     318           62 :      itypat=pawrhoij(iatom)%itypat
     319           62 :      mesh_size=pawtab(itypat)%mesh_size
     320           62 :      lm_size=pawtab(itypat)%lcut_size**2
     321           62 :      lm_size_eff=min(lm_size,pawang%ylm_size)
     322          186 :      ABI_MALLOC(ff,(mesh_size))
     323              : 
     324              : !    Compute on-site densities
     325          248 :      ABI_MALLOC(rho1 ,(mesh_size,lm_size,1))
     326          186 :      ABI_MALLOC(trho1,(mesh_size,lm_size,1))
     327          310 :      ABI_MALLOC(nhat1,(mesh_size,lm_size,usexcnhat))
     328          186 :      ABI_MALLOC(lmselect,(lm_size))
     329          124 :      ABI_MALLOC(lmselect_tmp,(lm_size))
     330          620 :      lmselect_tmp(:)=.true.
     331           62 :      cplex=1 ; nspden1=1 ; opt_compch=0 ; opt_dens=1-usexcnhat ; pawprtvol=0
     332              :      call pawdensities(compch_dum,cplex,iatom_tot,lmselect_tmp,lmselect,&
     333              : &     lm_size,nhat1,nspden1,nzlmopt,opt_compch,opt_dens,-1,0,pawang,pawprtvol,&
     334           62 : &     pawrad(itypat),pawrhoij(iatom),pawtab(itypat),rho1,trho1)
     335           62 :      ABI_FREE(lmselect_tmp)
     336              : 
     337              : !    Allocation of temporary densities
     338          124 :      ABI_MALLOC(rhotot,(mesh_size))
     339          186 :      ABI_MALLOC(grhotot,(mesh_size,3))
     340          186 :      ABI_MALLOC(rhoxc,(mesh_size,lm_size))
     341          124 :      ABI_MALLOC(drhoxc,(mesh_size))
     342          124 :      ABI_MALLOC(dcorexc,(mesh_size))
     343              : 
     344              : !    Loop over AE and PS contributions
     345          186 :      do iloop=1,2
     346          124 :        if (iloop==1) then
     347       776240 :          rhoxc(:,:)=rho1(:,:,1)
     348           62 :          corexc => pawtab(itypat)%coredens(:)
     349              :          usenhat=0
     350              :        else
     351       776240 :          rhoxc(:,:)=trho1(:,:,1)
     352           62 :          corexc => pawtab(itypat)%tcoredens(:,1)
     353              :          usenhat=usexcnhat
     354              :        end if
     355              : 
     356              : !      Add  hat density if needed
     357           62 :        if (usenhat>0) rhoxc(:,:)=rhoxc(:,:)+nhat1(:,:,1)
     358              : 
     359              : !      Need derivative of core density
     360              :        if (usecore==1) then
     361          124 :          call nderiv_gen(dcorexc,corexc,pawrad(itypat))
     362              :        end if
     363              : 
     364              : !      Do loop on the angular part (theta,phi)
     365        19530 :        do ipts=1,npts
     366              : 
     367              : !        Compute the density and its gradient for this (theta,phi)
     368    107630016 :          rhotot(:)=zero ; grhotot(:,:)=zero
     369       193440 :          do ilm=1,lm_size_eff
     370       193440 :            if (lmselect(ilm)) then
     371              :              !Density
     372     65966784 :              rhotot(:)=rhotot(:)+rhoxc(:,ilm)*pawang%ylmr(ilm,ipts)
     373              :              !Gradient
     374     65966784 :              ff(1:mesh_size)=rhoxc(1:mesh_size,ilm)
     375        47424 :              call nderiv_gen(drhoxc,ff,pawrad(itypat))
     376     65919360 :              ff(2:mesh_size)=ff(2:mesh_size)/pawrad(itypat)%rad(2:mesh_size)
     377        47424 :              call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
     378       189696 :              do ii=1,3
     379              :                grhotot(1:mesh_size,ii)=grhotot(1:mesh_size,ii) &
     380              : &                +drhoxc(1:mesh_size)*pawang%ylmr(ilm,ipts)*pawang%anginit(ii,ipts) &
     381    197947776 : &                +ff(1:mesh_size)*dylmdr(ii,ipts,ilm)
     382              :              end do
     383              :            end if
     384              :          end do
     385              :          if (usecore==1) then
     386     26907504 :            rhotot(:)=rhotot(:)+corexc(:)
     387        77376 :            do ii=1,3
     388     80741856 :              grhotot(:,ii)=grhotot(:,ii)+dcorexc(:)*pawang%anginit(ii,ipts)
     389              :            end do
     390              :          end if
     391              : 
     392              : !        Make the density positive everywhere (but do not care about gradients)
     393        19344 :          nspden1=1 ; iwarn=1
     394        19344 :          call mkdenpos(iwarn,mesh_size,nspden1,0,rhotot,xc_denpos)
     395              : 
     396              : !        Accumulate integral of |Grad_rho|/Rho
     397     26907504 :          do ir=1,mesh_size
     398     26888160 :            rho=rhotot(ir)
     399     26907504 :            if (abs(rho)>tol10) then
     400     26888160 :              grho2=grhotot(ir,1)**2+grhotot(ir,2)**2+grhotot(ir,3)**2
     401     26888160 :              ff(ir)=sqrt(grho2)/rho
     402              :            end if
     403              :          end do
     404     26907504 :          ff(1:mesh_size)=ff(1:mesh_size)*pawrad(itypat)%rad(1:mesh_size)**2
     405        19344 :          call simp_gen(sumg,ff,pawrad(itypat))
     406              : 
     407              : !        Add the contribution of the atom
     408        19468 :          if (iloop==1) then
     409         9672 :            grho_over_rho_paw=grho_over_rho_paw+sumg*four_pi*pawang%angwgth(ipts)
     410              :          else
     411         9672 :            grho_over_rho_paw=grho_over_rho_paw-sumg*four_pi*pawang%angwgth(ipts)
     412              :          end if
     413              : 
     414              :        end do ! npts
     415              :      end do ! AE/PS loop
     416              : 
     417              : !    Deallocate temporary memory
     418           62 :      ABI_FREE(ff)
     419           62 :      ABI_FREE(rho1)
     420           62 :      ABI_FREE(trho1)
     421           62 :      ABI_FREE(nhat1)
     422           62 :      ABI_FREE(lmselect)
     423           62 :      ABI_FREE(rhotot)
     424           62 :      ABI_FREE(grhotot)
     425           62 :      ABI_FREE(rhoxc)
     426           62 :      ABI_FREE(drhoxc)
     427           93 :      ABI_FREE(dcorexc)
     428              : 
     429              :    end do ! atom sites
     430              : 
     431              : !  Deallocate temporary memory
     432           31 :    ABI_FREE(dylmdr)
     433              : 
     434              : !  Reduction in case of parallelism
     435           31 :    if (paral_atom) then
     436            0 :      call xmpi_sum(grho_over_rho_paw,my_comm_atom,ierr)
     437              :    end if
     438              : 
     439              : !  Destroy atom table used for parallelism
     440           31 :    call free_my_atmtab(my_atmtab,my_atmtab_allocated)
     441              : 
     442              :  end if ! PAW and compute_paw
     443              : 
     444              : ! ========================================================================
     445              : ! ============ Assemble contributions and update mBJ C value  ============
     446              : 
     447           71 :  grho_over_rho = (grho_over_rho_pw + grho_over_rho_paw)/ucvol
     448           71 :  xc_c_tb09 = -0.012_dp + 1.023_dp * sqrt(grho_over_rho)
     449           71 :  if (xc_c_tb09<one) xc_c_tb09=one
     450           71 :  if (xc_c_tb09>two) xc_c_tb09=two
     451              : 
     452           71 :  write(msg,'(2a,f10.5,a)') ch10, &
     453          142 : &  ' In the mGGA functional TB09, c (computed value) = ',xc_c_tb09,ch10
     454           71 :  call wrtout(std_out,msg,'COLL')
     455              : 
     456           71 :  if (present(xc_funcs)) then
     457            0 :    call libxc_functionals_set_c_tb09(xc_c_tb09,xc_functionals=xc_funcs)
     458              :  else
     459           71 :    call libxc_functionals_set_c_tb09(xc_c_tb09)
     460              :  end if
     461              : 
     462              : !call timab(xx,2,tsec)
     463              : 
     464       110324 : end subroutine xc_tb09_update_c
     465              : !!***
     466              : 
     467              : !----------------------------------------------------------------------
     468              : 
     469              : END MODULE m_xc_tb09
     470              : !!***
        

Generated by: LCOV version 2.3-1