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 lda exchange hole in a truncated coulomb potential.
      10             : !>        Can be used as longrange correction for truncated hfx calculations
      11             : !> \par History
      12             : !>      Manuel Guidon (12.2008)  : created
      13             : !> \author Manuel Guidon (06.2008)
      14             : ! **************************************************************************************************
      15             : 
      16             : MODULE xc_xlda_hole_t_c_lr
      17             : 
      18             :    USE input_section_types,             ONLY: section_vals_type,&
      19             :                                               section_vals_val_get
      20             :    USE kinds,                           ONLY: dp
      21             :    USE mathconstants,                   ONLY: pi
      22             :    USE mathlib,                         ONLY: expint
      23             :    USE xc_derivative_desc,              ONLY: deriv_rho,&
      24             :                                               deriv_rhoa,&
      25             :                                               deriv_rhob
      26             :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      27             :                                               xc_dset_get_derivative
      28             :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      29             :                                               xc_derivative_type
      30             :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      31             :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      32             :                                               xc_rho_set_type
      33             : #include "../base/base_uses.f90"
      34             : 
      35             :    IMPLICIT NONE
      36             : 
      37             :    PRIVATE
      38             : 
      39             : ! *** Global parameters ***
      40             : 
      41             :    PUBLIC :: xlda_hole_t_c_lr_lda_eval, xlda_hole_t_c_lr_lda_info, &
      42             :              xlda_hole_t_c_lr_lsd_eval, xlda_hole_t_c_lr_lsd_info, &
      43             :              xlda_hole_t_c_lr_lda_calc_0
      44             : 
      45             :    REAL(KIND=dp), PARAMETER :: A = 1.0161144_dp, &
      46             :                                B = -0.37170836_dp, &
      47             :                                C = -0.077215461_dp, &
      48             :                                D = 0.57786348_dp, &
      49             :                                E = -0.051955731_dp
      50             : 
      51             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xlda_hole_t_c_lr'
      52             : 
      53             : CONTAINS
      54             : 
      55             : ! **************************************************************************************************
      56             : !> \brief returns various information on the functional
      57             : !> \param reference string with the reference of the actual functional
      58             : !> \param shortform string with the shortform of the functional name
      59             : !> \param needs the components needed by this functional are set to
      60             : !>        true (does not set the unneeded components to false)
      61             : !> \param max_deriv controls the number of derivatives
      62             : !> \par History
      63             : !>        12.2008 created [mguidon]
      64             : !> \author mguidon
      65             : ! **************************************************************************************************
      66          27 :    SUBROUTINE xlda_hole_t_c_lr_lda_info(reference, shortform, needs, max_deriv)
      67             :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      68             :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      69             :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      70             : 
      71          27 :       IF (PRESENT(reference)) THEN
      72           1 :          reference = "{LDA version}"
      73             :       END IF
      74          27 :       IF (PRESENT(shortform)) THEN
      75           1 :          shortform = "{LDA}"
      76             :       END IF
      77          27 :       IF (PRESENT(needs)) THEN
      78          26 :          needs%rho = .TRUE.
      79             :       END IF
      80          27 :       IF (PRESENT(max_deriv)) max_deriv = 1
      81             : 
      82          27 :    END SUBROUTINE xlda_hole_t_c_lr_lda_info
      83             : 
      84             : ! **************************************************************************************************
      85             : !> \brief returns various information on the functional
      86             : !> \param reference string with the reference of the actual functional
      87             : !> \param shortform string with the shortform of the functional name
      88             : !> \param needs the components needed by this functional are set to
      89             : !>        true (does not set the unneeded components to false)
      90             : !> \param max_deriv controls the number of derivatives
      91             : !> \par History
      92             : !>        12.2008 created [mguidon]
      93             : !> \author mguidon
      94             : ! **************************************************************************************************
      95          35 :    SUBROUTINE xlda_hole_t_c_lr_lsd_info(reference, shortform, needs, max_deriv)
      96             :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      97             :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      98             :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      99             : 
     100          35 :       IF (PRESENT(reference)) THEN
     101           1 :          reference = "{LSD version}"
     102             :       END IF
     103          35 :       IF (PRESENT(shortform)) THEN
     104           1 :          shortform = "{LSD}"
     105             :       END IF
     106          35 :       IF (PRESENT(needs)) THEN
     107          34 :          needs%rho_spin = .TRUE.
     108             :       END IF
     109          35 :       IF (PRESENT(max_deriv)) max_deriv = 1
     110             : 
     111          35 :    END SUBROUTINE xlda_hole_t_c_lr_lsd_info
     112             : 
     113             : ! **************************************************************************************************
     114             : !> \brief evaluates the truncated lda exchange hole
     115             : !> \param rho_set the density where you want to evaluate the functional
     116             : !> \param deriv_set place where to store the functional derivatives (they are
     117             : !>        added to the derivatives)
     118             : !> \param order degree of the derivative that should be evaluated,
     119             : !>        if positive all the derivatives up to the given degree are evaluated,
     120             : !>        if negative only the given degree is calculated
     121             : !> \param params input parameters (scaling, cutoff_radius)
     122             : !> \par History
     123             : !>      12.2008 created [Manuel Guidon]
     124             : !> \author Manuel Guidon
     125             : ! **************************************************************************************************
     126          88 :    SUBROUTINE xlda_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params)
     127             : 
     128             :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     129             :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     130             :       INTEGER, INTENT(IN)                                :: order
     131             :       TYPE(section_vals_type), POINTER                   :: params
     132             : 
     133             :       CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lda_eval'
     134             : 
     135             :       INTEGER                                            :: handle, npoints
     136             :       INTEGER, DIMENSION(2, 3)                           :: bo
     137             :       REAL(kind=dp)                                      :: epsilon_rho, R, sx
     138             :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     139          22 :          POINTER                                         :: dummy, e_0, e_rho, rho
     140             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     141             : 
     142          22 :       CALL timeset(routineN, handle)
     143             : 
     144          22 :       CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
     145          22 :       CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R)
     146             : 
     147             :       CALL xc_rho_set_get(rho_set, rho=rho, &
     148          22 :                           local_bounds=bo, rho_cutoff=epsilon_rho)
     149          22 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     150             : 
     151          22 :       dummy => rho
     152             : 
     153          22 :       e_0 => dummy
     154          22 :       e_rho => dummy
     155             : 
     156          22 :       IF (order >= 0) THEN
     157             :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     158          22 :                                          allocate_deriv=.TRUE.)
     159          22 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     160             :       END IF
     161          22 :       IF (order >= 1 .OR. order == -1) THEN
     162             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     163          22 :                                          allocate_deriv=.TRUE.)
     164          22 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     165             :       END IF
     166          22 :       IF (order > 1 .OR. order < -1) THEN
     167           0 :          CPABORT("derivatives bigger than 1 not implemented")
     168             :       END IF
     169             : 
     170          22 :       IF (R == 0.0_dp) THEN
     171           0 :          CPABORT("Cutoff_Radius 0.0 not implemented")
     172             :       END IF
     173             :       CALL xlda_hole_t_c_lr_lda_calc(npoints, order, rho=rho, &
     174             :                                      e_0=e_0, e_rho=e_rho, &
     175             :                                      epsilon_rho=epsilon_rho, &
     176          22 :                                      sx=sx, R=R)
     177             : 
     178          22 :       CALL timestop(handle)
     179             : 
     180          22 :    END SUBROUTINE xlda_hole_t_c_lr_lda_eval
     181             : 
     182             : ! **************************************************************************************************
     183             : !> \brief Call low level routine
     184             : !> \param npoints ...
     185             : !> \param order ...
     186             : !> \param rho ...
     187             : !> \param e_0 ...
     188             : !> \param e_rho ...
     189             : !> \param epsilon_rho ...
     190             : !> \param sx ...
     191             : !> \param R ...
     192             : !> \par History
     193             : !>      12.2008 created [Manuel Guidon]
     194             : !> \author Manuel Guidon
     195             : ! **************************************************************************************************
     196          22 :    SUBROUTINE xlda_hole_t_c_lr_lda_calc(npoints, order, rho, e_0, e_rho, &
     197             :                                         epsilon_rho, sx, R)
     198             : 
     199             :       INTEGER, INTENT(in)                                :: npoints, order
     200             :       REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
     201             :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho, sx, R
     202             : 
     203             :       INTEGER                                            :: ip
     204             :       REAL(dp)                                           :: my_rho
     205             : 
     206             : !$OMP     PARALLEL DO DEFAULT(NONE) &
     207             : !$OMP                 SHARED(npoints, rho, epsilon_rho, order, e_0, e_rho) &
     208             : !$OMP                 SHARED(sx, r) &
     209          22 : !$OMP                 PRIVATE(ip, my_rho)
     210             : 
     211             :       DO ip = 1, npoints
     212             :          my_rho = MAX(rho(ip), 0.0_dp)
     213             :          IF (my_rho > epsilon_rho) THEN
     214             :             CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_0(ip), e_rho(ip), &
     215             :                                              sx, R)
     216             :          END IF
     217             :       END DO
     218             : 
     219             : !$OMP     END PARALLEL DO
     220             : 
     221          22 :    END SUBROUTINE xlda_hole_t_c_lr_lda_calc
     222             : 
     223             : ! **************************************************************************************************
     224             : !> \brief low level routine
     225             : !> \param order ...
     226             : !> \param rho ...
     227             : !> \param e_0 ...
     228             : !> \param e_rho ...
     229             : !> \param sx ...
     230             : !> \param R ...
     231             : !> \par History
     232             : !>      12.2008 created [Manuel Guidon]
     233             : !> \author Manuel Guidon
     234             : ! **************************************************************************************************
     235     9091946 :    SUBROUTINE xlda_hole_t_c_lr_lda_calc_0(order, rho, e_0, e_rho, &
     236             :                                           sx, R)
     237             :       INTEGER, INTENT(IN)                                :: order
     238             :       REAL(KIND=dp), INTENT(IN)                          :: rho
     239             :       REAL(kind=dp), INTENT(INOUT)                       :: e_0, e_rho
     240             :       REAL(KIND=dp), INTENT(IN)                          :: sx, R
     241             : 
     242             :       REAL(KIND=dp)                                      :: t1, t12, t14, t15, t19, t2, t22, t23, &
     243             :                                                             t24, t25, t3, t32, t33, t36, t4, t41, &
     244             :                                                             t46, t5, t6, t62, t64, t67, t68, t7, &
     245             :                                                             t82, t86, t9, t91, t95
     246             : 
     247     9091946 :       IF (order >= 0) THEN
     248     9091946 :          t1 = rho**2
     249     9091946 :          t2 = t1*pi
     250     9091946 :          t3 = 3**(0.1e1_dp/0.3e1_dp)
     251     9091946 :          t4 = pi**2
     252     9091946 :          t5 = t4*rho
     253     9091946 :          t6 = t5**(0.1e1_dp/0.3e1_dp)
     254     9091946 :          t7 = t6**2
     255     9091946 :          t9 = t3/t7
     256     9091946 :          t12 = LOG(R*t3*t6)
     257     9091946 :          t14 = R**2
     258     9091946 :          t15 = t14**2
     259     9091946 :          t19 = 0.1e1_dp/D
     260     9091946 :          t22 = t3**2
     261     9091946 :          t23 = t22*t7
     262     9091946 :          t24 = D*t14*t23
     263     9091946 :          t25 = EXP(-t24)
     264     9091946 :          t32 = 9 + 4*A*t14*t23
     265     9091946 :          t33 = LOG(t32)
     266     9091946 :          t36 = D**2
     267     9091946 :          t41 = expint(1, t24)
     268     9091946 :          t46 = 0.1e1_dp/t36
     269     9091946 :          t62 = LOG(0.2e1_dp)
     270     9091946 :          t64 = LOG(A)
     271             :          t67 = A*t12 + 0.3e1_dp/0.2e1_dp*E*t15*t3*t6*t5*t19*t25 &
     272             :                - A*t33/0.2e1_dp + E/t36/D*t25 + A*t41/0.2e1_dp + E*t14 &
     273             :                *t22*t7*t46*t25 + B*t19*t25/0.2e1_dp + C*t46*t25/0.2e1_dp &
     274             :                + C*t14*t22*t7*t19*t25/0.2e1_dp + A*t62 + A*t64 &
     275     9091946 :                /0.2e1_dp
     276     9091946 :          t68 = t9*t67
     277     9091946 :          e_0 = e_0 + (0.2e1_dp/0.3e1_dp*t2*t68)*sx
     278             :       END IF
     279     9091946 :       IF (order >= 1 .OR. order == -1) THEN
     280     4940448 :          t82 = A/rho
     281     4940448 :          t86 = t4**2
     282     4940448 :          t91 = A**2
     283     4940448 :          t95 = 0.1e1_dp/t6*t4
     284             :          e_rho = e_rho + (0.4e1_dp/0.3e1_dp*rho*pi*t68 - 0.4e1_dp/0.9e1_dp*t1*t4*pi &
     285             :                           *t3/t7/t5*t67 + 0.2e1_dp/0.3e1_dp*t2*t9*(t82/0.3e1_dp - &
     286             :                                                                    0.3e1_dp*E*t15*t14*t86*rho*t25 - 0.4e1_dp/0.3e1_dp*t91*t14 &
     287             :                                                                    *t22*t95/t32 - t82*t25/0.3e1_dp - B*t14*t22*t95*t25 &
     288     4940448 :                                                                    /0.3e1_dp - C*t15*t3*t6*t4*t25))*sx
     289             :       END IF
     290             : 
     291     9091946 :    END SUBROUTINE xlda_hole_t_c_lr_lda_calc_0
     292             : 
     293             : ! **************************************************************************************************
     294             : !> \brief evaluates the truncated lsd exchange hole. Calls the lda routine and
     295             : !>        applies spin scaling relation
     296             : !> \param rho_set the density where you want to evaluate the functional
     297             : !> \param deriv_set place where to store the functional derivatives (they are
     298             : !>        added to the derivatives)
     299             : !> \param order degree of the derivative that should be evaluated,
     300             : !>        if positive all the derivatives up to the given degree are evaluated,
     301             : !>        if negative only the given degree is calculated
     302             : !> \param params input parameters (scaling, cutoff_radius)
     303             : !> \par History
     304             : !>      12.2008 created [Manuel Guidon]
     305             : !> \author Manuel Guidon
     306             : ! **************************************************************************************************
     307         120 :    SUBROUTINE xlda_hole_t_c_lr_lsd_eval(rho_set, deriv_set, order, params)
     308             : 
     309             :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     310             :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     311             :       INTEGER, INTENT(IN)                                :: order
     312             :       TYPE(section_vals_type), POINTER                   :: params
     313             : 
     314             :       CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lsd_eval'
     315             : 
     316             :       INTEGER                                            :: handle, npoints
     317             :       INTEGER, DIMENSION(2, 3)                           :: bo
     318             :       REAL(kind=dp)                                      :: epsilon_rho, R, sx
     319             :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     320          30 :          POINTER                                         :: dummy, e_0, e_rhoa, e_rhob, rhoa, rhob
     321             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     322             : 
     323          30 :       CALL timeset(routineN, handle)
     324             : 
     325          30 :       CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
     326          30 :       CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R)
     327             : 
     328             :       CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, &
     329          30 :                           local_bounds=bo, rho_cutoff=epsilon_rho)
     330          30 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     331             : 
     332          30 :       dummy => rhoa
     333             : 
     334          30 :       e_0 => dummy
     335          30 :       e_rhoa => dummy
     336          30 :       e_rhob => dummy
     337             : 
     338          30 :       IF (order >= 0) THEN
     339             :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     340          30 :                                          allocate_deriv=.TRUE.)
     341          30 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     342             :       END IF
     343          30 :       IF (order >= 1 .OR. order == -1) THEN
     344             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
     345          16 :                                          allocate_deriv=.TRUE.)
     346          16 :          CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
     347             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
     348          16 :                                          allocate_deriv=.TRUE.)
     349          16 :          CALL xc_derivative_get(deriv, deriv_data=e_rhob)
     350             :       END IF
     351          30 :       IF (order > 1 .OR. order < -1) THEN
     352           0 :          CPABORT("derivatives bigger than 2 not implemented")
     353             :       END IF
     354          30 :       IF (R == 0.0_dp) THEN
     355           0 :          CPABORT("Cutoff_Radius 0.0 not implemented")
     356             :       END IF
     357             : 
     358             : !$OMP     PARALLEL DEFAULT(NONE) &
     359             : !$OMP              SHARED(npoints, order, rhoa, e_0, e_rhoa, epsilon_rho) &
     360          30 : !$OMP              SHARED(sx, r,rhob, e_rhob)
     361             : 
     362             :       CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhoa, &
     363             :                                      e_0=e_0, e_rho=e_rhoa, &
     364             :                                      epsilon_rho=epsilon_rho, &
     365             :                                      sx=sx, R=R)
     366             : 
     367             :       CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhob, &
     368             :                                      e_0=e_0, e_rho=e_rhob, &
     369             :                                      epsilon_rho=epsilon_rho, &
     370             :                                      sx=sx, R=R)
     371             : !$OMP     END PARALLEL
     372             : 
     373          30 :       CALL timestop(handle)
     374             : 
     375          30 :    END SUBROUTINE xlda_hole_t_c_lr_lsd_eval
     376             : 
     377             : ! **************************************************************************************************
     378             : !> \brief low level routine
     379             : !> \param npoints ...
     380             : !> \param order ...
     381             : !> \param rho ...
     382             : !> \param e_0 ...
     383             : !> \param e_rho ...
     384             : !> \param epsilon_rho ...
     385             : !> \param sx ...
     386             : !> \param R ...
     387             : !> \par History
     388             : !>      12.2008 created [Manuel Guidon]
     389             : !> \author Manuel Guidon
     390             : ! **************************************************************************************************
     391          60 :    SUBROUTINE xlda_hole_t_c_lr_lsd_calc(npoints, order, rho, e_0, e_rho, &
     392             :                                         epsilon_rho, sx, R)
     393             : 
     394             :       INTEGER, INTENT(in)                                :: npoints, order
     395             :       REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
     396             :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho, sx, R
     397             : 
     398             :       INTEGER                                            :: ip
     399             :       REAL(dp)                                           :: e_tmp, my_rho
     400             : 
     401             : !$OMP     DO
     402             : 
     403             :       DO ip = 1, npoints
     404     2733750 :          my_rho = 2.0_dp*MAX(rho(ip), 0.0_dp)
     405     2733750 :          IF (my_rho > epsilon_rho) THEN
     406     2733750 :             e_tmp = 0.0_dp
     407             :             CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_tmp, e_rho(ip), &
     408     2733750 :                                              sx, R)
     409     2733750 :             e_0(ip) = e_0(ip) + 0.5_dp*e_tmp
     410             :          END IF
     411             :       END DO
     412             : 
     413             : !$OMP     END DO
     414             : 
     415          60 :    END SUBROUTINE xlda_hole_t_c_lr_lsd_calc
     416             : END MODULE xc_xlda_hole_t_c_lr
     417             :