LCOV - code coverage report
Current view: top level - src/xc - xc_xbeef.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:f7c1873) Lines: 68 148 45.9 %
Date: 2024-03-26 06:29:31 Functions: 3 6 50.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief calculates the Exchange contribution in the BEEF-vdW functional
      10             : !         (Norskov Group)
      11             : !> \par History
      12             : !>      02.2014 created based on xc_xbecke88.F [rkoitz]
      13             : !> \author rkoitz
      14             : ! **************************************************************************************************
      15             : MODULE xc_xbeef
      16             : 
      17             :    USE bibliography,                    ONLY: Wellendorff2012,&
      18             :                                               cite_reference
      19             :    USE cp_array_utils,                  ONLY: cp_3d_r_cp_type
      20             :    USE input_section_types,             ONLY: section_vals_type,&
      21             :                                               section_vals_val_get
      22             :    USE kinds,                           ONLY: dp
      23             :    USE mathconstants,                   ONLY: pi
      24             :    USE xc_derivative_desc,              ONLY: deriv_norm_drho,&
      25             :                                               deriv_norm_drhoa,&
      26             :                                               deriv_norm_drhob,&
      27             :                                               deriv_rho,&
      28             :                                               deriv_rhoa,&
      29             :                                               deriv_rhob
      30             :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      31             :                                               xc_dset_get_derivative
      32             :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      33             :                                               xc_derivative_type
      34             :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      35             :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      36             :                                               xc_rho_set_type
      37             : #include "../base/base_uses.f90"
      38             : 
      39             :    IMPLICIT NONE
      40             :    PRIVATE
      41             : 
      42             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      43             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbeef'
      44             : 
      45             :    PUBLIC :: xbeef_lda_info, xbeef_lsd_info, xbeef_lda_eval, xbeef_lsd_eval
      46             : 
      47             :    !for beef functional pulled out of GPAW source, Feb 2014
      48             :    REAL(kind=dp), PARAMETER :: a(0:29) = (/1.516501714304992365356_dp, 0.441353209874497942611_dp, -0.091821352411060291887_dp, &
      49             :                                            -0.023527543314744041314_dp, 0.034188284548603550816_dp, 0.002411870075717384172_dp, &
      50             :                                            -0.014163813515916020766_dp, 0.000697589558149178113_dp, 0.009859205136982565273_dp, &
      51             :                                            -0.006737855050935187551_dp, -0.001573330824338589097_dp, 0.005036146253345903309_dp, &
      52             :                                            -0.002569472452841069059_dp, -0.000987495397608761146_dp, 0.002033722894696920677_dp, &
      53             :                                            -0.000801871884834044583_dp, -0.000668807872347525591_dp, 0.001030936331268264214_dp, &
      54             :                                            -0.000367383865990214423_dp, -0.000421363539352619543_dp, 0.000576160799160517858_dp, &
      55             :                                            -0.000083465037349510408_dp, -0.000445844758523195788_dp, 0.000460129009232047457_dp, &
      56             :                                            -0.000005231775398304339_dp, -0.000423957047149510404_dp, 0.000375019067938866537_dp, &
      57             :                                            0.000021149381251344578_dp, -0.000190491156503997170_dp, 0.000073843624209823442_dp/)
      58             : CONTAINS
      59             : 
      60             : ! **************************************************************************************************
      61             : !> \brief return various information on the functional
      62             : !> \param reference string with the reference of the actual functional
      63             : !> \param shortform string with the shortform of the functional name
      64             : !> \param needs the components needed by this functional are set to
      65             : !>        true (does not set the unneeded components to false)
      66             : !> \param max_deriv ...
      67             : !> \par History
      68             : !>      02.2014 created
      69             : !> \author rkoitz
      70             : ! **************************************************************************************************
      71             : 
      72          23 :    SUBROUTINE xbeef_lda_info(reference, shortform, needs, max_deriv)
      73             :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      74             :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      75             :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      76             : 
      77          23 :       IF (PRESENT(reference)) THEN
      78           1 :          reference = "Wellendorff, J. et al., Phys. Rev. B 85, 235149 (2012) {LDA}"
      79             :       END IF
      80          23 :       IF (PRESENT(shortform)) THEN
      81           1 :          shortform = "Exchange Contribution to BEEF-vdW Functional (Wellendorff, 2012) {LDA}"
      82             :       END IF
      83          23 :       IF (PRESENT(needs)) THEN
      84          22 :          needs%rho = .TRUE.
      85          22 :          needs%rho_1_3 = .TRUE.
      86          22 :          needs%norm_drho = .TRUE.
      87             :       END IF
      88          23 :       IF (PRESENT(max_deriv)) max_deriv = 1
      89             : 
      90          23 :    END SUBROUTINE xbeef_lda_info
      91             : 
      92             : ! **************************************************************************************************
      93             : !> \brief return various information on the functional
      94             : !> \param reference string with the reference of the actual functional
      95             : !> \param shortform string with the shortform of the functional name
      96             : !> \param needs the components needed by this functional are set to
      97             : !>        true (does not set the unneeded components to false)
      98             : !> \param max_deriv ...
      99             : !> \par History
     100             : !>      02.2014 created
     101             : !> \author rkoitz
     102             : ! **************************************************************************************************
     103           0 :    SUBROUTINE xbeef_lsd_info(reference, shortform, needs, max_deriv)
     104             :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
     105             :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
     106             :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
     107             : 
     108           0 :       IF (PRESENT(reference)) THEN
     109           0 :          reference = "Wellendorff, J. et al., Phys. Rev. B 85, 235149 (2012) {LSD}"
     110             :       END IF
     111           0 :       IF (PRESENT(shortform)) THEN
     112           0 :          shortform = "Exchange Contribution to BEEF-vdW Functional (Wellendorff, 2012) {LSD}"
     113             :       END IF
     114           0 :       IF (PRESENT(needs)) THEN
     115           0 :          needs%rho_spin = .TRUE.
     116           0 :          needs%rho_spin_1_3 = .TRUE.
     117           0 :          needs%norm_drho_spin = .TRUE.
     118             :       END IF
     119           0 :       IF (PRESENT(max_deriv)) max_deriv = 1
     120             : 
     121           0 :    END SUBROUTINE xbeef_lsd_info
     122             : 
     123             : ! **************************************************************************************************
     124             : !> \brief evaluates the beef exchange functional for lda
     125             : !> \param rho_set the density where you want to evaluate the functional
     126             : !> \param deriv_set place where to store the functional derivatives (they are
     127             : !>        added to the derivatives)
     128             : !> \param grad_deriv degree of the derivative that should be evaluated,
     129             : !>        if positive all the derivatives up to the given degree are evaluated,
     130             : !>        if negative only the given degree is calculated
     131             : !> \param xbeef_params input parameters (scaling)
     132             : !> \par History
     133             : !>      02.2014 created
     134             : !> \author rkoitz
     135             : ! **************************************************************************************************
     136          54 :    SUBROUTINE xbeef_lda_eval(rho_set, deriv_set, grad_deriv, xbeef_params)
     137             :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     138             :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     139             :       INTEGER, INTENT(in)                                :: grad_deriv
     140             :       TYPE(section_vals_type), POINTER                   :: xbeef_params
     141             : 
     142             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xbeef_lda_eval'
     143             : 
     144             :       INTEGER                                            :: handle, npoints
     145             :       INTEGER, DIMENSION(2, 3)                           :: bo
     146             :       REAL(kind=dp)                                      :: epsilon_rho, sx
     147             :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     148          18 :          POINTER                                         :: dummy, e_0, e_ndrho, e_rho, norm_drho, &
     149          18 :                                                             rho, rho_1_3
     150             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     151             : 
     152          18 :       CALL timeset(routineN, handle)
     153             : 
     154          18 :       CALL section_vals_val_get(xbeef_params, "scale_x", r_val=sx)
     155             : 
     156          18 :       CALL cite_reference(Wellendorff2012)
     157             : 
     158             :       CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, rho=rho, &
     159          18 :                           norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
     160          18 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     161             : 
     162          18 :       dummy => rho
     163             : 
     164          18 :       e_0 => dummy
     165          18 :       e_rho => dummy
     166          18 :       e_ndrho => dummy
     167             : 
     168          18 :       IF (grad_deriv >= 0) THEN
     169             :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     170          18 :                                          allocate_deriv=.TRUE.)
     171          18 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     172             :       END IF
     173          18 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     174             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     175          18 :                                          allocate_deriv=.TRUE.)
     176          18 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     177             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     178          18 :                                          allocate_deriv=.TRUE.)
     179          18 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     180             :       END IF
     181          18 :       IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
     182           0 :          CPABORT("derivatives greater than 1 not implemented")
     183             :       END IF
     184             : 
     185             : !$OMP     PARALLEL DEFAULT(NONE) &
     186             : !$OMP              SHARED(rho, rho_1_3, norm_drho, e_0, e_rho) &
     187             : !$OMP              SHARED(e_ndrho) &
     188             : !$OMP              SHARED( grad_deriv, npoints) &
     189          18 : !$OMP              SHARED(epsilon_rho,sx)
     190             :       CALL xbeef_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
     191             :                           e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
     192             :                           grad_deriv=grad_deriv, &
     193             :                           npoints=npoints, epsilon_rho=epsilon_rho, sx=sx)
     194             : !$OMP     END PARALLEL
     195             : 
     196          18 :       CALL timestop(handle)
     197          18 :    END SUBROUTINE xbeef_lda_eval
     198             : 
     199             : ! **************************************************************************************************
     200             : !> \brief evaluates the beef exchange functional for lda
     201             : !> \param rho the density where you want to evaluate the functional
     202             : !> \param rho_1_3 ...
     203             : !> \param norm_drho ...
     204             : !> \param e_0 ...
     205             : !> \param e_rho ...
     206             : !> \param e_ndrho ...
     207             : !> \param grad_deriv degree of the derivative that should be evaluated,
     208             : !>        if positive all the derivatives up to the given degree are evaluated,
     209             : !>        if negative only the given degree is calculated
     210             : !> \param npoints ...
     211             : !> \param epsilon_rho ...
     212             : !> \param sx scaling-parameter for exchange
     213             : !> \par History
     214             : !>      02.2014 created based on xc_xbecke88
     215             : !> \author rkoitz
     216             : ! **************************************************************************************************
     217          18 :    SUBROUTINE xbeef_lda_calc(rho, rho_1_3, norm_drho, &
     218          18 :                              e_0, e_rho, e_ndrho, &
     219             :                              grad_deriv, npoints, epsilon_rho, sx)
     220             :       INTEGER, INTENT(in)                                :: npoints, grad_deriv
     221             :       REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho, e_rho, e_0
     222             :       REAL(kind=dp), DIMENSION(1:npoints), INTENT(in)    :: norm_drho, rho_1_3, rho
     223             :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho, sx
     224             : 
     225             :       INTEGER, PARAMETER                                 :: m = 30
     226             : 
     227             :       INTEGER                                            :: i, ii
     228             :       REAL(kind=dp)                                      :: ds_ndrho, ds_rho, dt, e_ueg, e_ueg_drho, &
     229             :                                                             epsilon_rho43, k_lda, kf, my_rho, &
     230             :                                                             my_rho_1_3, s, s2, t, t3
     231             :       REAL(kind=dp), DIMENSION(0:29)                     :: de_leg, e_leg
     232             : 
     233             : !energies, first and second derivatives from legendre pol
     234             : 
     235          18 :       kf = (3.0_dp*pi**2)**(1.0_dp/3.0_dp) !only constants, without n^1/3
     236          18 :       k_lda = -(2.0_dp/(2.0_dp**(4._dp/3._dp)))*(3.0_dp/2.0_dp)*(3.0_dp/(4.0_dp*pi))**(1.0_dp/3.0_dp)
     237             :       !exchange energy density of the uniform electron gas, non spin-polarized
     238             : 
     239          18 :       epsilon_rho43 = epsilon_rho**(4._dp/3._dp)
     240             : 
     241             : !$OMP     DO
     242             :       DO ii = 1, npoints
     243             : 
     244     1417176 :          my_rho = rho(ii)
     245             : 
     246     1417176 :          IF (my_rho > epsilon_rho) THEN
     247     1417176 :             my_rho_1_3 = rho_1_3(ii)
     248             : 
     249     1417176 :             e_ueg = k_lda*my_rho*my_rho_1_3
     250     1417176 :             e_ueg_drho = (4.0_dp/3.0_dp)*k_lda*my_rho_1_3
     251             : 
     252     1417176 :             t3 = my_rho_1_3*my_rho*2*kf !reduced gradient, denominator
     253             : 
     254     1417176 :             s = norm_drho(ii)/MAX(t3, epsilon_rho43) !reduced gradient finally
     255     1417176 :             s2 = s**2
     256     1417176 :             t = 2.0_dp*s2/(4.0_dp + s2) - 1.0_dp
     257             : 
     258     1417176 :             IF (grad_deriv >= 0) THEN !asking for pure e evaluation or also derivatives
     259     1417176 :                e_leg(0) = 1 !first legendre pol
     260     1417176 :                e_leg(1) = t !second legendre pol
     261             :             END IF
     262             : 
     263     1417176 :             IF ((grad_deriv >= 1) .OR. (grad_deriv == -1)) THEN !asking for first derivative or higher
     264     1417176 :                de_leg(0) = 0
     265     1417176 :                de_leg(1) = 1
     266     1417176 :                dt = 4.0_dp*s/(4.0_dp + s2) - 4.0_dp*s*s2/(4.0_dp + s2)**2
     267     1417176 :                ds_rho = -(4.0_dp*s)/(3.0_dp*MAX(my_rho, epsilon_rho))
     268     1417176 :                ds_ndrho = 1.0_dp/(MAX(t3, epsilon_rho43))
     269             :             END IF
     270             : 
     271    41098104 :             DO i = 2, m - 1 !LEGENDRE PART
     272    39680928 :                e_leg(i) = 2.*(t)*e_leg(i - 1) - e_leg(i - 2) - ((t)*e_leg(i - 1) - e_leg(i - 2))/(REAL(i, KIND=dp))
     273             :                !taken from quantum espresso beef library.
     274             : 
     275    41098104 :                IF (ABS(grad_deriv) >= 1) THEN !first derivative
     276             :                   !the zero-derivatives need to be available for the first deriv.
     277    39680928 :                   de_leg(i) = e_leg(i - 1)*i + de_leg(i - 1)*(t)
     278             :                END IF
     279             :             END DO
     280             : 
     281             :             !NO DERIVATIVE
     282     1417176 :             IF (grad_deriv >= 0) THEN
     283             :                !add the scaled legendre linear combination to e_0
     284    43932456 :                e_0(ii) = e_0(ii) + SUM(e_leg*a)*e_ueg*sx
     285             :             END IF
     286             : 
     287             :             !FIRST DERIVATIVE
     288     1417176 :             IF ((grad_deriv >= 1) .OR. (grad_deriv == -1)) THEN !asking for first derivative or higher
     289    86447736 :                e_rho(ii) = e_rho(ii) + (SUM(e_leg*a)*e_ueg_drho + SUM(de_leg*a)*dt*ds_rho*e_ueg)*sx
     290    43932456 :                e_ndrho(ii) = e_ndrho(ii) + (SUM(de_leg*a)*dt*ds_ndrho*e_ueg)*sx
     291             :             END IF
     292             : 
     293             :          END IF
     294             :       END DO
     295             : 
     296             : !$OMP    END DO
     297             : 
     298          18 :    END SUBROUTINE xbeef_lda_calc
     299             : 
     300             : ! **************************************************************************************************
     301             : !> \brief evaluates the beef 88 exchange functional for lsd
     302             : !> \param rho_set ...
     303             : !> \param deriv_set ...
     304             : !> \param grad_deriv ...
     305             : !> \param xbeef_params ...
     306             : !> \par History
     307             : !>         2/2014 rkoitz [created based on Becke 88]
     308             : !> \author rkoitz
     309             : ! **************************************************************************************************
     310           0 :    SUBROUTINE xbeef_lsd_eval(rho_set, deriv_set, grad_deriv, xbeef_params)
     311             :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     312             :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     313             :       INTEGER, INTENT(in)                                :: grad_deriv
     314             :       TYPE(section_vals_type), POINTER                   :: xbeef_params
     315             : 
     316             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xbeef_lsd_eval'
     317             : 
     318             :       INTEGER                                            :: handle, i, ispin, npoints
     319             :       INTEGER, DIMENSION(2, 3)                           :: bo
     320             :       REAL(kind=dp)                                      :: epsilon_rho, sx
     321             :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     322           0 :          POINTER                                         :: dummy, e_0
     323           0 :       TYPE(cp_3d_r_cp_type), DIMENSION(2)                :: e_ndrho, e_rho, norm_drho, rho, rho_1_3
     324             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     325             : 
     326           0 :       CALL timeset(routineN, handle)
     327             : 
     328           0 :       CALL cite_reference(Wellendorff2012)
     329             : 
     330           0 :       NULLIFY (deriv)
     331           0 :       DO i = 1, 2
     332           0 :          NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
     333             :       END DO
     334             : 
     335           0 :       CALL section_vals_val_get(xbeef_params, "scale_x", r_val=sx)
     336             :       CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
     337             :                           rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
     338             :                           rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
     339             :                           norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
     340           0 :                           local_bounds=bo)
     341           0 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     342             : 
     343           0 :       dummy => rho(1)%array
     344             : 
     345           0 :       e_0 => dummy
     346           0 :       DO i = 1, 2
     347           0 :          e_rho(i)%array => dummy
     348           0 :          e_ndrho(i)%array => dummy
     349             :       END DO
     350             : 
     351           0 :       IF (grad_deriv >= 0) THEN
     352             :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     353           0 :                                          allocate_deriv=.TRUE.)
     354           0 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     355             :       END IF
     356           0 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     357             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
     358           0 :                                          allocate_deriv=.TRUE.)
     359           0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
     360             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
     361           0 :                                          allocate_deriv=.TRUE.)
     362           0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
     363             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
     364           0 :                                          allocate_deriv=.TRUE.)
     365           0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
     366             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
     367           0 :                                          allocate_deriv=.TRUE.)
     368           0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
     369             :       END IF
     370           0 :       IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
     371           0 :          CPABORT("derivatives greater than 1 not implemented")
     372             :       END IF
     373             : 
     374           0 :       DO ispin = 1, 2
     375             : 
     376             : !$OMP        PARALLEL DEFAULT(NONE) &
     377             : !$OMP                 SHARED(rho, ispin, rho_1_3, norm_drho, e_0) &
     378             : !$OMP                 SHARED(e_rho, e_ndrho) &
     379             : !$OMP                 SHARED(grad_deriv, npoints) &
     380           0 : !$OMP                 SHARED(epsilon_rho, sx)
     381             : 
     382             :          CALL xbeef_lsd_calc( &
     383             :             rho_spin=rho(ispin)%array, &
     384             :             rho_1_3_spin=rho_1_3(ispin)%array, &
     385             :             norm_drho_spin=norm_drho(ispin)%array, &
     386             :             e_0=e_0, e_rho_spin=e_rho(ispin)%array, &
     387             :             e_ndrho_spin=e_ndrho(ispin)%array, &
     388             :             grad_deriv=grad_deriv, npoints=npoints, &
     389             :             epsilon_rho=epsilon_rho, sx=sx)
     390             : 
     391             : !$OMP        END PARALLEL
     392             : 
     393             :       END DO
     394             : 
     395           0 :       CALL timestop(handle)
     396             : 
     397           0 :    END SUBROUTINE xbeef_lsd_eval
     398             : ! **************************************************************************************************
     399             : !> \brief low level calculation of the beef exchange functional for lsd
     400             : !> \param rho_spin alpha or beta spin density
     401             : !> \param rho_1_3_spin rho_spin**(1./3.)
     402             : !> \param norm_drho_spin || grad rho_spin ||
     403             : !> \param e_0 adds to it the local value of the functional
     404             : !> \param e_rho_spin e_*_spin: derivative of the functional wrt. to the variables
     405             : !>        named where the * is. Everything wrt. to the spin of the arguments.
     406             : !> \param e_ndrho_spin ...
     407             : !> \param grad_deriv ...
     408             : !> \param npoints ...
     409             : !> \param epsilon_rho ...
     410             : !> \param sx scaling-parameter for exchange
     411             : !> \par History
     412             : !>      02.2014 created based on Becke88
     413             : !> \author rkoitz
     414             : ! **************************************************************************************************
     415           0 :    SUBROUTINE xbeef_lsd_calc(rho_spin, rho_1_3_spin, norm_drho_spin, e_0, &
     416             :                              e_rho_spin, e_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx)
     417             :       REAL(kind=dp), DIMENSION(*), INTENT(in)            :: rho_spin, rho_1_3_spin, norm_drho_spin
     418             :       REAL(kind=dp), DIMENSION(*), INTENT(inout)         :: e_0, e_rho_spin, e_ndrho_spin
     419             :       INTEGER, INTENT(in)                                :: grad_deriv, npoints
     420             :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho, sx
     421             : 
     422             :       INTEGER, PARAMETER                                 :: m = 30
     423             : 
     424             :       INTEGER                                            :: i, ii
     425             :       REAL(kind=dp)                                      :: ds_ndrho, ds_rho, dt, e_ueg, e_ueg_drho, &
     426             :                                                             epsilon_rho43, k_lsd, kf, &
     427             :                                                             my_epsilon_rho, my_rho, my_rho_1_3, s, &
     428             :                                                             s2, t, t3
     429             :       REAL(kind=dp), DIMENSION(0:29)                     :: de_leg, e_leg
     430             : 
     431             : !energies and first derivatives from legendre pol
     432             : 
     433           0 :       kf = (3.0_dp*pi**2)**(1.0_dp/3.0_dp) !only constants, without n^1/3
     434           0 :       k_lsd = (3.0_dp/2.0_dp)*(3.0_dp/(4.0_dp*pi))**(1.0_dp/3.0_dp)
     435             :       !exchange energy density of the uniform electron gas, spin-polarized
     436             : 
     437           0 :       my_epsilon_rho = 0.5_dp*epsilon_rho
     438           0 :       epsilon_rho43 = my_epsilon_rho**(4._dp/3._dp)
     439             : 
     440           0 : !$OMP     DO
     441             :       DO ii = 1, npoints
     442           0 :          my_rho = rho_spin(ii)
     443             : 
     444           0 :          IF (my_rho > epsilon_rho) THEN
     445           0 :             my_rho_1_3 = rho_1_3_spin(ii)
     446             : 
     447           0 :             e_ueg = k_lsd*my_rho*my_rho_1_3
     448           0 :             e_ueg_drho = (4.0_dp/3.0_dp)*k_lsd*my_rho_1_3
     449             : 
     450           0 :             t3 = my_rho_1_3*my_rho*2*kf !reduced gradient, denominator
     451             : 
     452           0 :             s = norm_drho_spin(ii)/MAX(t3, epsilon_rho43) !reduced gradient finally
     453           0 :             s2 = s**2
     454           0 :             t = 2.0_dp*s**2/(4.0_dp + s**2) - 1.0_dp
     455             : 
     456           0 :             IF (grad_deriv >= 0) THEN !asking for pure e evaluation or also derivatives
     457           0 :                e_leg(0) = 1 !first legendre pol
     458           0 :                e_leg(1) = t !second legendre pol
     459             :             END IF
     460             : 
     461           0 :             IF ((grad_deriv >= 1) .OR. (grad_deriv == -1)) THEN !asking for first derivative or higher
     462           0 :                de_leg(0) = 0
     463           0 :                de_leg(1) = 1
     464           0 :                dt = 4.0_dp*s/(4.0_dp + s2) - 4.0_dp*s*s2/(4.0_dp + s2)**2
     465           0 :                ds_rho = -(4.0_dp*s)/(3.0_dp*MAX(my_rho, epsilon_rho))
     466           0 :                ds_ndrho = 1.0_dp/(MAX(t3, epsilon_rho43))
     467             :             END IF
     468             : 
     469           0 :             DO i = 2, m - 1 !LEGENDRE PART
     470           0 :                e_leg(i) = 2.*(t)*e_leg(i - 1) - e_leg(i - 2) - ((t)*e_leg(i - 1) - e_leg(i - 2))/(REAL(i, KIND=dp))
     471             :                !taken from quantum espresso beef library.
     472             : 
     473           0 :                IF (ABS(grad_deriv) >= 1) THEN !first derivative
     474             :                   !the zero-derivatives need to be available for the first deriv.
     475           0 :                   de_leg(i) = e_leg(i - 1)*i + de_leg(i - 1)*(t)
     476             :                END IF
     477             : 
     478             :             END DO
     479             : 
     480             :             !NO DERIVATIVE
     481           0 :             IF (grad_deriv >= 0) THEN
     482             :                !add the scaled legendre linear combination to e_0
     483           0 :                e_0(ii) = e_0(ii) + SUM(e_leg*a)*e_ueg*sx
     484             :             END IF
     485             : 
     486             :             !FIRST DERIVATIVE
     487           0 :             IF ((grad_deriv >= 1) .OR. (grad_deriv == -1)) THEN !asking for first derivative or higher
     488           0 :                e_rho_spin(ii) = e_rho_spin(ii) + (SUM(e_leg*a)*e_ueg_drho + SUM(de_leg*a)*dt*ds_rho*e_ueg)*sx
     489           0 :                e_ndrho_spin(ii) = e_ndrho_spin(ii) + (SUM(de_leg*a)*dt*ds_ndrho*e_ueg)*sx
     490             :             END IF
     491             :          END IF
     492             :       END DO
     493             : !$OMP   END DO
     494             : 
     495           0 :    END SUBROUTINE xbeef_lsd_calc
     496             : 
     497             : END MODULE xc_xbeef
     498             : 

Generated by: LCOV version 1.15