LCOV - code coverage report
Current view: top level - types - types_xcpot_libxc.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 151 203 74.4 %
Date: 2019-09-08 04:53:50 Functions: 19 22 86.4 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 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             : 
       7             : !>This module contains the xcpot-type providing an interface to libxc
       8             : MODULE m_types_xcpot_libxc
       9             : #ifdef CPP_LIBXC
      10             :    USE xc_f03_lib_m
      11             : #endif
      12             :    USE m_types_xcpot
      13             :    USE m_judft
      14             :    use m_types_misc
      15             :    IMPLICIT NONE
      16             : 
      17             : 
      18             : #ifdef CPP_LIBXC
      19             :    PRIVATE :: write_xc_info
      20             : #endif
      21             : 
      22             :    TYPE,EXTENDS(t_xcpot):: t_xcpot_libxc 
      23             : #ifdef CPP_LIBXC    
      24             :       TYPE(xc_f03_func_t)      :: vxc_func_x, vxc_func_c
      25             :       TYPE(xc_f03_func_t)      :: exc_func_x, exc_func_c
      26             : #endif
      27             :       INTEGER                  :: func_vxc_id_c, func_vxc_id_x !> functionals to be used for potential & density convergence
      28             :       INTEGER                  :: func_exc_id_c, func_exc_id_x !> functionals to be used in exc- & totale-calculations
      29             :       INTEGER                  :: jspins
      30             :    CONTAINS
      31             :       PROCEDURE        :: vx_is_LDA     => xcpot_vx_is_LDA
      32             :       PROCEDURE        :: vx_is_GGA     => xcpot_vx_is_GGA
      33             :       PROCEDURE        :: vx_is_MetaGGA => xcpot_vx_is_MetaGGA
      34             :       PROCEDURE        :: vx_is_XC      => xcpot_vx_is_XC ! exchange and correlation in one functional
      35             :       
      36             :       PROCEDURE        :: vc_is_LDA => xcpot_vc_is_LDA
      37             :       PROCEDURE        :: vc_is_GGA => xcpot_vc_is_GGA
      38             : 
      39             :       PROCEDURE        :: exc_is_LDA     => xcpot_exc_is_LDA
      40             :       PROCEDURE        :: exc_is_gga     => xcpot_exc_is_gga
      41             :       PROCEDURE        :: exc_is_MetaGGA => xcpot_exc_is_MetaGGA
      42             :       PROCEDURE        :: ex_is_XC       => xcpot_ex_is_XC ! exchange and correlation in one functional
      43             : 
      44             :       PROCEDURE        :: is_hybrid           => xcpot_is_hybrid 
      45             :       PROCEDURE        :: get_exchange_weight => xcpot_get_exchange_weight
      46             :       PROCEDURE        :: get_vxc             => xcpot_get_vxc
      47             :       PROCEDURE        :: get_exc             => xcpot_get_exc
      48             :       PROCEDURE,NOPASS :: alloc_gradients     => xcpot_alloc_gradients
      49             :       !Not             overloeaded...
      50             :       PROCEDURE        :: init                => xcpot_init
      51             :    END TYPE t_xcpot_libxc
      52             :    PUBLIC t_xcpot_libxc
      53             : CONTAINS
      54             : 
      55           4 :    SUBROUTINE xcpot_init(xcpot,jspins,vxc_id_x,vxc_id_c,exc_id_x,exc_id_c)
      56             :       USE m_judft
      57             :       IMPLICIT NONE
      58             :    CLASS(t_xcpot_libxc),INTENT(OUT)    :: xcpot
      59             :       INTEGER, INTENT(IN)                 :: jspins
      60             :       INTEGER, INTENT(IN)                 :: vxc_id_x, vxc_id_c ! potential functional
      61             :       INTEGER, INTENT(IN)                 :: exc_id_x, exc_id_c ! energy functionals
      62             :       LOGICAL                             :: same_functionals   ! are vxc and exc equal
      63             :       INTEGER                             :: errors(4)
      64             : 
      65             : #ifdef CPP_LIBXC
      66             :       INTEGER :: err
      67           4 :       xcpot%jspins = jspins
      68             : 
      69           4 :       xcpot%func_vxc_id_x = vxc_id_x
      70           4 :       xcpot%func_vxc_id_c = vxc_id_c
      71             :       
      72           4 :       xcpot%func_exc_id_x = exc_id_x
      73           4 :       xcpot%func_exc_id_c = exc_id_c
      74             : 
      75           4 :       if(xcpot%func_vxc_id_x == 0 .or. xcpot%func_exc_id_x == 0 ) then
      76             :          CALL judft_error("LibXC exchange- and correlation-function indicies need to be set"&
      77             :             ,hint='Try this: ' // ACHAR(10) //&
      78             :             '<xcFunctional name="libxc" relativisticCorrections="F">' // ACHAR(10) //& 
      79             :             '  <libXC  exchange="1" correlation="1" /> ' // ACHAR(10) //& 
      80           0 :             '</xcFunctional> ')
      81             :       endif 
      82             : 
      83           4 :       IF (jspins==1) THEN
      84             :          ! potential functionals      
      85           2 :          CALL xc_f03_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_UNPOLARIZED, err=errors(1))
      86           2 :          IF (xcpot%func_vxc_id_c>0) CALL xc_f03_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, &
      87           2 :                                                                  XC_UNPOLARIZED, err=errors(2))
      88             : 
      89             :          ! energy functionals
      90           2 :          CALL xc_f03_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_UNPOLARIZED, err=errors(3))
      91           2 :          IF (xcpot%func_exc_id_c>0) CALL xc_f03_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, &
      92           2 :                                                                   XC_UNPOLARIZED, err=errors(4))
      93             :          
      94             :       ELSE
      95             :          ! potential functionals  
      96           2 :          CALL xc_f03_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_POLARIZED, err=errors(1))
      97           2 :          IF (xcpot%func_vxc_id_c>0) CALL xc_f03_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, &
      98           2 :                                                                   XC_POLARIZED, err=errors(2))
      99             : 
     100             :          !energy functionals
     101           2 :          CALL xc_f03_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_POLARIZED, err=errors(3))
     102           2 :          IF (xcpot%func_exc_id_c>0) CALL xc_f03_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, &
     103           2 :                                                                   XC_POLARIZED, err=errors(4))
     104             :       END IF
     105             : 
     106           4 :       IF(errors(1) /= 0) call juDFT_error("Exchange potential functional not in LibXC")
     107           4 :       IF(errors(2) /= 0) call juDFT_error("Correlation potential functional not in LibXC")
     108           4 :       IF(errors(3) /= 0) call juDFT_error("Exchange energy functional not in LibXC")
     109           4 :       IF(errors(4) /= 0) call juDFT_error("Correlation energy functional not in LibXC")
     110             : 
     111             :       !check if any potental is a MetaGGA
     112           8 :       IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_c))) THEN
     113           0 :          call juDFT_error("vxc_x: MetaGGA is not implemented for correclation potentials")
     114             :       ENDIF
     115             : 
     116           4 :       CALL write_xc_info(xcpot%vxc_func_x)
     117             : 
     118           4 :       IF (xcpot%func_vxc_id_c>0) THEN
     119           4 :          CALL write_xc_info(xcpot%vxc_func_c)
     120             :       ELSE
     121           0 :          WRITE(*,*) "No Correlation functional"
     122             :       END IF
     123             : 
     124             :       same_functionals =     (xcpot%func_vxc_id_x == xcpot%func_exc_id_x) &
     125           4 :                        .and. (xcpot%func_vxc_id_c == xcpot%func_exc_id_c)
     126           4 :       IF(.not. same_functionals) THEN
     127           0 :          CALL write_xc_info(xcpot%exc_func_x)
     128           0 :          IF (xcpot%func_exc_id_c>0) THEN
     129           0 :             CALL write_xc_info(xcpot%exc_func_c)
     130             :          ELSE
     131           0 :             WRITE(*,*) "No Correlation functional for TotalE"
     132             :          ENDIF
     133             :       ELSE
     134           4 :          write (*,*) "Using same functional for VXC and EXC"
     135             :       END IF
     136             : 
     137           4 :       xcpot%kinED%set = .False.
     138             : #else
     139             :       CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
     140             :          hint="Please recompile FLEUR with libxc support")
     141             : #endif
     142           4 :    END SUBROUTINE xcpot_init
     143             : 
     144             :    ! LDA
     145           2 :    LOGICAL FUNCTION xcpot_vx_is_LDA(xcpot)
     146             :       IMPLICIT NONE
     147             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     148             : #ifdef CPP_LIBXC
     149             :       TYPE(xc_f03_func_info_t)        :: xc_info
     150             : 
     151           2 :       xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
     152           2 :       xcpot_vx_is_LDA =  XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info)
     153             : #else
     154             :       xcpot_vx_is_LDA = .false.
     155             : #endif
     156           2 :    END FUNCTION xcpot_vx_is_LDA
     157             :    
     158           1 :    LOGICAL FUNCTION xcpot_vc_is_LDA(xcpot)
     159             :       IMPLICIT NONE
     160             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     161             : #ifdef CPP_LIBXC
     162             :       TYPE(xc_f03_func_info_t)        :: xc_info
     163             : 
     164           1 :       xc_info = xc_f03_func_get_info(xcpot%vxc_func_c)
     165           1 :       xcpot_vc_is_LDA =  XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info)
     166             : #else
     167             :       xcpot_vc_is_LDA = .false.
     168             : #endif
     169           1 :    END FUNCTION xcpot_vc_is_LDA
     170             :    
     171           3 :    LOGICAL FUNCTION xcpot_exc_is_LDA(xcpot)
     172             :       IMPLICIT NONE
     173             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     174             : #ifdef CPP_LIBXC
     175             :       TYPE(xc_f03_func_info_t)        :: xc_info
     176             : 
     177           3 :       xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
     178           3 :       xcpot_exc_is_LDA = (XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info))
     179             : #else
     180             :       xcpot_exc_is_LDA = .false.
     181             : #endif
     182           3 :    END FUNCTION xcpot_exc_is_LDA
     183             : 
     184             :    ! GGA
     185      118300 :    LOGICAL FUNCTION xcpot_vc_is_gga(xcpot)
     186             :       IMPLICIT NONE
     187             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     188             : #ifdef CPP_LIBXC
     189             :       TYPE(xc_f03_func_info_t)        :: xc_info
     190             : 
     191      118300 :       xc_info = xc_f03_func_get_info(xcpot%vxc_func_c)
     192      118300 :       xcpot_vc_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
     193             : #else
     194             :       xcpot_vc_is_gga=.false.
     195             : #endif
     196      118300 :    END FUNCTION xcpot_vc_is_gga
     197             : 
     198           1 :    LOGICAL FUNCTION xcpot_vx_is_gga(xcpot)
     199             :       IMPLICIT NONE
     200             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     201             : #ifdef CPP_LIBXC
     202             :       TYPE(xc_f03_func_info_t)        :: xc_info
     203             : 
     204           1 :       xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
     205           1 :       xcpot_vx_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
     206             : #else
     207             :       xcpot_vx_is_gga=.false.
     208             : #endif
     209           1 :    END FUNCTION xcpot_vx_is_gga
     210             : 
     211        2465 :    LOGICAL FUNCTION xcpot_vx_is_MetaGGA(xcpot)
     212             :       IMPLICIT NONE
     213             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     214             : #ifdef CPP_LIBXC
     215             :       TYPE(xc_f03_func_info_t)        :: xc_info
     216             : 
     217        2465 :       xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
     218        2465 :       xcpot_vx_is_MetaGGA =  ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f03_func_info_get_family(xc_info))
     219             : #else
     220             :       xcpot_vx_is_MetaGGA=.false.
     221             : #endif
     222        2465 :    END FUNCTION xcpot_vx_is_MetaGGA
     223             :    
     224         110 :    LOGICAL FUNCTION xcpot_vx_is_XC(xcpot)
     225             :       IMPLICIT NONE
     226             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     227             : #ifdef CPP_LIBXC
     228             :       TYPE(xc_f03_func_info_t)        :: xc_info
     229             : 
     230         110 :       xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
     231         110 :       xcpot_vx_is_XC =  xc_f03_func_info_get_kind(xc_info)  == XC_EXCHANGE_CORRELATION
     232             : #else
     233             :       xcpot_vx_is_XC=.false.
     234             : #endif
     235         110 :    END FUNCTION xcpot_vx_is_XC
     236             :    
     237          63 :    LOGICAL FUNCTION xcpot_ex_is_XC(xcpot)
     238             :       IMPLICIT NONE
     239             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     240             : #ifdef CPP_LIBXC
     241             :       TYPE(xc_f03_func_info_t)        :: xc_info
     242             : 
     243          63 :       xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
     244          63 :       xcpot_ex_is_XC =  xc_f03_func_info_get_kind(xc_info)  == XC_EXCHANGE_CORRELATION
     245             : #else
     246             :       xcpot_ex_is_XC=.false.
     247             : #endif
     248          63 :    END FUNCTION xcpot_ex_is_XC
     249             : 
     250          63 :    LOGICAL FUNCTION xcpot_exc_is_gga(xcpot)
     251             :       IMPLICIT NONE
     252             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     253             : #ifdef CPP_LIBXC
     254             :       TYPE(xc_f03_func_info_t)        :: xc_info
     255             : 
     256          63 :       xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
     257          63 :       xcpot_exc_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
     258             : #else
     259             :       xcpot_exc_is_gga=.false.
     260             : #endif
     261          63 :    END FUNCTION xcpot_exc_is_gga
     262             : 
     263         184 :    LOGICAL FUNCTION xcpot_exc_is_MetaGGA(xcpot)
     264             :       IMPLICIT NONE
     265             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     266             : #ifdef CPP_LIBXC    
     267             :       TYPE(xc_f03_func_info_t)        :: xc_info
     268             : 
     269         184 :       xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
     270         184 :       xcpot_exc_is_MetaGGA=ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f03_func_info_get_family(xc_info))
     271             : #else
     272             :       xcpot_exc_is_MetaGGA =  .False.
     273             : #endif
     274         184 :    END FUNCTION xcpot_exc_is_MetaGGA
     275             : 
     276          50 :    LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
     277             :       IMPLICIT NONE
     278             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     279             : #ifdef CPP_LIBXC
     280             :       TYPE(xc_f03_func_info_t)        :: xc_info
     281             : 
     282          50 :       xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
     283          50 :       xcpot_is_hybrid=ANY([XC_FAMILY_HYB_MGGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
     284             : #else
     285             :       xcpot_is_hybrid = .False.
     286             : #endif
     287          50 :    END FUNCTION xcpot_is_hybrid
     288             : 
     289           0 :    FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
     290             :       USE m_judft
     291             :       IMPLICIT NONE
     292             :    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
     293             : 
     294             :       REAL:: a_ex
     295             : #ifdef CPP_LIBXC    
     296           0 :       a_ex=xc_f03_hyb_exx_coef(xcpot%vxc_func_x)
     297             : #endif
     298           0 :    END FUNCTION xcpot_get_exchange_weight
     299             : 
     300             :    !***********************************************************************
     301         220 :    SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad, kinED_KS)
     302             :       USE, INTRINSIC :: IEEE_ARITHMETIC
     303             :       IMPLICIT NONE
     304             :    CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
     305             :       INTEGER, INTENT (IN)     :: jspins
     306             :       REAL,INTENT (IN)         :: rh(:,:)   !Dimensions here
     307             :       REAL, INTENT (OUT)       :: vx (:,:)  !points,spin
     308             :       REAL, INTENT (OUT  )     :: vxc(:,:)  !
     309             :       REAL, INTENT(IN),OPTIONAL:: kinED_KS(:,:)
     310             :       ! optional arguments for GGA
     311             :       TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
     312             : #ifdef CPP_LIBXC    
     313         110 :       REAL,ALLOCATABLE  :: vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:), &
     314         110 :                            tmp_vsig(:,:), tmp_vlapl(:,:), tmp_vtau(:,:), &
     315         110 :                            kinED_libxc(:,:)
     316             :       integer           :: idx
     317             :       !libxc uses the spin as a first index, hence we have to transpose....
     318         110 :       ALLOCATE(vxc_tmp(SIZE(vxc,2),SIZE(vxc,1)));vxc_tmp=0.0
     319         110 :       ALLOCATE(vx_tmp(SIZE(vx,2),SIZE(vx,1)));vx_tmp=0.0
     320             : 
     321         110 :       IF (xcpot%vc_is_GGA()) THEN
     322          85 :          IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
     323          85 :          ALLOCATE(vsigma,mold=grad%vsigma)
     324             :          !where(abs(grad%sigma)<1E-9) grad%sigma=1E-9
     325          85 :          CALL xc_f03_gga_vxc(xcpot%vxc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vx_tmp,vsigma)
     326             : 
     327             :          ! check if vx has both exchange and correlation
     328          85 :          IF (.not. xcpot%vx_is_XC()) THEN
     329          85 :             CALL xc_f03_gga_vxc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vxc_tmp,grad%vsigma)
     330          85 :             grad%vsigma=grad%vsigma+vsigma
     331          85 :             vxc_tmp=vxc_tmp+vx_tmp
     332             :          ENDIF
     333             :       ELSE  !LDA potentials
     334          25 :          if(xcpot%vx_is_MetaGGA() .and. present(kinED_KS)) then
     335           0 :             kinED_libXC = transpose(kinED_KS + 0.25 * grad%laplace)
     336             : 
     337           0 :             allocate(tmp_vsig,  mold=vx_tmp)
     338           0 :             allocate(tmp_vlapl, mold=vx_tmp)
     339           0 :             allocate(tmp_vtau,  mold=vx_tmp)
     340             :             
     341             :             call xc_f03_mgga_vxc(xcpot%vxc_func_x, size(rh,1), transpose(rh), &
     342             :                                  grad%sigma, transpose(grad%laplace), kinED_libxc,&
     343           0 :                                  vx_tmp, tmp_vsig, tmp_vlapl, tmp_vtau)
     344             : 
     345           0 :             idx = find_first_normal(vx_tmp)+1
     346           0 :             vx_tmp(:,:idx) = 0.0
     347             : 
     348           0 :             CALL xc_f03_lda_vxc(initial_lda_func(jspins), idx, TRANSPOSE(rh(:idx,:)), vx_tmp(:,:idx))
     349             :          else
     350          25 :             CALL xc_f03_lda_vxc(initial_lda_func(jspins), SIZE(rh,1), TRANSPOSE(rh), vx_tmp)
     351             :          endif
     352          25 :          IF (.not. xcpot%vx_is_XC()) THEN
     353          25 :             CALL xc_f03_lda_vxc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh), vxc_tmp)
     354          25 :             vxc_tmp=vxc_tmp+vx_tmp
     355             :          ENDIF
     356             :       ENDIF
     357         110 :       vx=TRANSPOSE(vx_tmp) 
     358         305 :       vxc=TRANSPOSE(vxc_tmp) 
     359             : #endif
     360         110 :    END SUBROUTINE xcpot_get_vxc
     361             : 
     362           0 :    FUNCTION find_first_normal(vec) result(idx)
     363         110 :       USE, INTRINSIC :: IEEE_ARITHMETIC
     364             :       implicit none
     365             :       real, intent(in)    :: vec(:,:)
     366             :       integer             :: idx
     367             : 
     368           0 :       idx = size(vec, dim=2)
     369           0 :       do while(all(.not. ieee_is_nan(vec(:,idx))))
     370           0 :          idx = idx - 1
     371           0 :          if(idx == 0) exit
     372             :       enddo
     373           0 :    END FUNCTION find_first_normal
     374             : 
     375             : 
     376          63 :    SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinED_KS, mt_call)
     377           0 :       use m_constants
     378             :       IMPLICIT NONE
     379             :    CLASS(t_xcpot_libxc),INTENT(IN)          :: xcpot
     380             :       INTEGER, INTENT (IN)                  :: jspins
     381             :       REAL,INTENT (IN)                      :: rh(:,:)  !points,spin
     382             :       REAL, INTENT (OUT)                    :: exc(:) !points
     383             :       ! optional arguments for GGA
     384             :       TYPE(t_gradients),OPTIONAL,INTENT(IN) :: grad
     385             :       LOGICAL, OPTIONAL, INTENT(IN)         :: mt_call    
     386             : 
     387             :       ! kinED from Kohn-Sham equations:
     388             :       ! tau = sum[phi_i(r)^dag nabla phi_i(r)]
     389             :       ! see eq (2) in https://doi.org/10.1063/1.1565316
     390             :       ! (-0.5 is applied below)
     391             :       REAL, INTENT(IN), OPTIONAL     :: kinED_KS(:,:)
     392             : 
     393             : #ifdef CPP_LIBXC
     394             :       TYPE(xc_f03_func_info_t)       :: xc_info
     395         126 :       REAL  :: excc(SIZE(exc))
     396             :       REAL  :: cut_ratio = 0.1
     397             :       INTEGER :: cut_idx
     398             :       LOGICAL :: is_mt
     399             : 
     400             :       ! tau = 0.5 * sum[|grad phi_i(r)|²]
     401             :       ! see eq (3) in https://doi.org/10.1063/1.1565316
     402          63 :       REAL, ALLOCATABLE              :: kinED_libXC(:,:), pkzb_ratio(:,:), pkzb_zaehler(:,:), pkzb_nenner(:,:)
     403             : 
     404          63 :       is_mt = merge(mt_call, .False., present(mt_call))
     405          63 :       IF (xcpot%exc_is_gga()) THEN
     406          60 :          IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a GGA potential without providing derivatives")
     407          60 :          CALL xc_f03_gga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
     408             : 
     409             :          !in case ex is both exchange and correlation
     410          60 :          IF (.not. xcpot%ex_is_XC()) THEN
     411          60 :             CALL xc_f03_gga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc)
     412     3944460 :             exc=exc+excc
     413             :          END IF
     414           3 :       ELSEIF(xcpot%exc_is_LDA()) THEN  !LDA potentials
     415           3 :          CALL xc_f03_lda_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc)
     416             :          
     417             :          !in case ex is both exchange and correlation
     418           3 :          IF (.not. xcpot%ex_is_XC()) THEN
     419           3 :             CALL xc_f03_lda_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), excc)
     420      706775 :             exc=exc+excc
     421             :          END IF
     422           0 :       ELSEIF(xcpot%exc_is_MetaGGA()) THEN
     423           0 :          IF(PRESENT(kinED_KS)) THEN 
     424             :             ! apply correction in  eq (4) in https://doi.org/10.1063/1.1565316
     425           0 :             kinED_libXC = transpose(kinED_KS + 0.25 * grad%laplace)
     426             : 
     427             :             !only cut core of muffin tin
     428           0 :             cut_idx = MERGE(NINT(size(rh,1) * cut_ratio), 0, is_mt)
     429             : 
     430           0 :             exc  = 0.0
     431           0 :             excc = 0.0
     432             :             call xc_f03_mgga_exc(xcpot%exc_func_x, SIZE(rh(cut_idx+1:,:),1),& 
     433             :                                                    TRANSPOSE(rh(cut_idx+1:,:)), &
     434             :                                                    grad%sigma(:,cut_idx+1:), &
     435             :                                                    transpose(grad%laplace(cut_idx+1:,:)), &
     436             :                                                    kinED_libXC(:,cut_idx+1:), &
     437           0 :                                                    exc(cut_idx+1:))
     438             : 
     439             :             call xc_f03_gga_exc(xcpot%vxc_func_x, SIZE(rh(:cut_idx,:),1), &
     440             :                                                   TRANSPOSE(rh(:cut_idx,:)), &
     441             :                                                   grad%sigma(:,:cut_idx), &
     442           0 :                                                   exc(:cut_idx))
     443             : 
     444           0 :             IF (xcpot%func_exc_id_c>0) THEN
     445             :                call xc_f03_mgga_exc(xcpot%exc_func_c, SIZE(rh(cut_idx+1:,:),1), &
     446             :                                                       TRANSPOSE(rh(cut_idx+1:,:)), &
     447             :                                                       grad%sigma(:,cut_idx+1:), &
     448             :                                                       transpose(grad%laplace(cut_idx+1:,:)), &
     449             :                                                       kinED_libXC(:,cut_idx+1:), &
     450           0 :                                                       excc(cut_idx+1:))
     451             :                
     452             :                call xc_f03_gga_exc(xcpot%vxc_func_c, SIZE(rh(:cut_idx,:),1), &
     453             :                                                       TRANSPOSE(rh(:cut_idx,:)), &
     454             :                                                       grad%sigma(:,:cut_idx), &
     455           0 :                                                       excc(:cut_idx))
     456           0 :                exc=exc+excc
     457             :             END IF
     458             :             
     459             :          ELSE ! first iteration is GGA
     460           0 :             IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a MetaGGA potential without providing derivatives")
     461           0 :             CALL xc_f03_gga_exc(xcpot%vxc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
     462           0 :             IF (xcpot%func_exc_id_c>0) THEN
     463           0 :                CALL xc_f03_gga_exc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc)
     464           0 :                exc=exc+excc
     465             :             END IF
     466             :          ENDIF
     467             : 
     468             :       ELSE
     469           0 :          call juDFT_error("exc is part of a known Family", calledby="xcpot_get_exc@libxc")
     470             :       ENDIF
     471             : 
     472             : #endif
     473         126 :    END SUBROUTINE xcpot_get_exc
     474             : 
     475         105 :    SUBROUTINE xcpot_alloc_gradients(ngrid,jspins,grad)
     476             :       INTEGER, INTENT (IN)         :: jspins,ngrid
     477             :       TYPE(t_gradients),INTENT(INOUT):: grad
     478             :       !For libxc we only need the sigma array...
     479         105 :       IF (ALLOCATED(grad%sigma)) DEALLOCATE(grad%sigma,grad%gr,grad%laplace,grad%vsigma)
     480         105 :       ALLOCATE(grad%sigma(MERGE(1,3,jspins==1),ngrid))
     481         105 :       ALLOCATE(grad%gr(3,ngrid,jspins))
     482         105 :       ALLOCATE(grad%laplace(ngrid,jspins))
     483         105 :       ALLOCATE(grad%vsigma(MERGE(1,3,jspins==1),ngrid))
     484             : 
     485         335 :       grad%gr      = 0.0
     486         335 :       grad%laplace = 0.0
     487     4063705 :       grad%sigma   = 0.0
     488     4063705 :       grad%vsigma  = 0.0
     489         105 :    END SUBROUTINE xcpot_alloc_gradients
     490             : 
     491             : #ifdef CPP_LIBXC  
     492           8 :    SUBROUTINE write_xc_info(xc_func, is_E_func)
     493             :       IMPLICIT NONE
     494             :       LOGICAL,INTENT(IN),OPTIONAL         :: is_E_func
     495             :       INTEGER                             :: i
     496             :       CHARACTER(len=120)                  :: kind, family
     497             :       LOGICAL                             :: is_energy_func
     498             : 
     499             :       TYPE(xc_f03_func_t),INTENT(IN)      :: xc_func
     500             :       TYPE(xc_f03_func_info_t)            :: xc_info
     501             : 
     502           8 :       xc_info = xc_f03_func_get_info(xc_func)
     503           8 :       is_energy_func = merge(is_E_func, .False., PRESENT(is_E_func))
     504             : 
     505           4 :       SELECT CASE(xc_f03_func_info_get_kind(xc_info))
     506             :       CASE (XC_EXCHANGE)
     507           4 :          WRITE(kind, '(a)') 'an exchange functional'
     508             :       CASE (XC_CORRELATION)
     509           4 :          WRITE(kind, '(a)') 'a correlation functional'
     510             :       CASE (XC_EXCHANGE_CORRELATION)
     511           0 :          WRITE(kind, '(a)') 'an exchange-correlation functional'
     512             :       CASE (XC_KINETIC)
     513           0 :          WRITE(kind, '(a)') 'a kinetic energy functional'
     514             :       CASE default
     515           8 :          WRITE(kind, '(a)') 'of unknown kind'
     516             :       END SELECT
     517           4 :       SELECT CASE (xc_f03_func_info_get_family(xc_info))
     518             :       CASE (XC_FAMILY_LDA);
     519           4 :          WRITE(family,'(a)') "LDA"
     520             :       CASE (XC_FAMILY_GGA);
     521           4 :          WRITE(family,'(a)') "GGA"
     522             :       CASE (XC_FAMILY_HYB_GGA);
     523           0 :          WRITE(family,'(a)') "Hybrid GGA"
     524             :       CASE (XC_FAMILY_MGGA);
     525           0 :          WRITE(family,'(a)') "MGGA"
     526             :       CASE (XC_FAMILY_HYB_MGGA);
     527           0 :          WRITE(family,'(a)') "Hybrid MGGA"
     528             :       CASE default;
     529           8 :          WRITE(family,'(a)') "unknown"
     530             :       END SELECT
     531             : 
     532           8 :       IF(.not. is_energy_func) THEN
     533             :          WRITE(*,'("The functional ''", a, "'' is ", a, ", it belongs to the ''", a, "'' family and is defined in the reference(s):")') &
     534           8 :             TRIM(xc_f03_func_info_get_name(xc_info)), TRIM(kind), TRIM(family)
     535             :       ELSE
     536             :          WRITE(*,'("The functional used for TotalE ''", a, "'' is ", a, ", it belongs to the ''", a, "'' family and is defined in the reference(s):")') &
     537           0 :             TRIM(xc_f03_func_info_get_name(xc_info)), TRIM(kind), TRIM(family)
     538             :       ENDIF
     539             : 
     540           8 :       i = 0
     541          36 :       DO WHILE(i >= 0)
     542          22 :          WRITE(*, '(a,i1,2a)') '[', i+1, '] ', TRIM(xc_f03_func_reference_get_ref(xc_f03_func_info_get_references(xc_info, i)))
     543             :       END DO
     544           8 :    END SUBROUTINE write_xc_info
     545             : 
     546           0 :    FUNCTION xc_get_family(xc_func) result(family)
     547             :       IMPLICIT NONE
     548             :       TYPE(xc_f03_func_t)  :: xc_func  
     549             :       integer              :: family
     550           4 :       family = xc_f03_func_info_get_family(xc_f03_func_get_info(xc_func))
     551           0 :    END FUNCTION xc_get_family
     552             : 
     553          25 :    function initial_lda_func(jspins) result(lda)
     554             :       use, intrinsic :: iso_c_binding
     555             :       implicit none
     556             :       INTEGER, intent(in) :: jspins
     557             :       TYPE(xc_f03_func_t) :: lda
     558             :       integer             :: ierr
     559             : 
     560          25 :       if(jspins == 1) then
     561          25 :          CALL xc_f03_func_init(lda, 1, XC_UNPOLARIZED, err=ierr)
     562             :       else
     563           0 :          CALL xc_f03_func_init(lda, 1, XC_POLARIZED, err=ierr)
     564             :       ENDIF
     565             : 
     566          25 :       if(ierr /= 0) call juDFT_error("can't find lda_x for init")
     567             : 
     568          50 :    end function initial_lda_func
     569             : 
     570             : #endif
     571             : 
     572           8 : END MODULE m_types_xcpot_libxc

Generated by: LCOV version 1.13