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 calculate optx
      10             : !> \note
      11             : !>      will need proper testing / review
      12             : !> \author Joost VandeVondele [03.2004]
      13             : ! **************************************************************************************************
      14             : MODULE xc_optx
      15             :    USE cp_array_utils,                  ONLY: cp_3d_r_cp_type
      16             :    USE input_section_types,             ONLY: section_vals_type,&
      17             :                                               section_vals_val_get
      18             :    USE kinds,                           ONLY: dp
      19             :    USE xc_derivative_desc,              ONLY: deriv_norm_drho,&
      20             :                                               deriv_norm_drhoa,&
      21             :                                               deriv_norm_drhob,&
      22             :                                               deriv_rho,&
      23             :                                               deriv_rhoa,&
      24             :                                               deriv_rhob
      25             :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      26             :                                               xc_dset_get_derivative
      27             :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      28             :                                               xc_derivative_type
      29             :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      30             :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      31             :                                               xc_rho_set_type
      32             : #include "../base/base_uses.f90"
      33             : 
      34             :    IMPLICIT NONE
      35             :    PRIVATE
      36             : 
      37             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_optx'
      38             : 
      39             :    PUBLIC :: optx_lda_info, optx_lda_eval, optx_lsd_info, optx_lsd_eval
      40             : CONTAINS
      41             : 
      42             : ! **************************************************************************************************
      43             : !> \brief info about the optx functional
      44             : !> \param reference string with the reference of the actual functional
      45             : !> \param shortform string with the shortform of the functional name
      46             : !> \param needs the components needed by this functional are set to
      47             : !>        true (does not set the unneeded components to false)
      48             : !> \param max_deriv implemented derivative of the xc functional
      49             : !> \author Joost
      50             : ! **************************************************************************************************
      51         435 :    SUBROUTINE optx_lda_info(reference, shortform, needs, max_deriv)
      52             :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      53             :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      54             :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      55             : 
      56         435 :       IF (PRESENT(reference)) THEN
      57           3 :          reference = "OPTX, Handy NC and Cohen AJ,  JCP 116, p. 5411 (2002) (LDA)"
      58             :       END IF
      59         435 :       IF (PRESENT(shortform)) THEN
      60           3 :          shortform = "OPTX exchange (LDA)"
      61             :       END IF
      62         435 :       IF (PRESENT(needs)) THEN
      63         432 :          needs%rho = .TRUE.
      64         432 :          needs%norm_drho = .TRUE.
      65             :       END IF
      66         435 :       IF (PRESENT(max_deriv)) max_deriv = 1
      67         435 :    END SUBROUTINE optx_lda_info
      68             : 
      69             : ! **************************************************************************************************
      70             : !> \brief info about the optx functional (LSD)
      71             : !> \param reference string with the reference of the actual functional
      72             : !> \param shortform string with the shortform of the functional name
      73             : !> \param needs the components needed by this functional are set to
      74             : !>        true (does not set the unneeded components to false)
      75             : !> \param max_deriv implemented derivative of the xc functional
      76             : !> \author Joost
      77             : ! **************************************************************************************************
      78         473 :    SUBROUTINE optx_lsd_info(reference, shortform, needs, max_deriv)
      79             :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      80             :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      81             :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      82             : 
      83         473 :       IF (PRESENT(reference)) THEN
      84           1 :          reference = "OPTX, Handy NC and Cohen AJ,  JCP 116, p. 5411 (2002), (LSD) "
      85             :       END IF
      86         473 :       IF (PRESENT(shortform)) THEN
      87           1 :          shortform = "OPTX exchange (LSD)"
      88             :       END IF
      89         473 :       IF (PRESENT(needs)) THEN
      90         472 :          needs%rho_spin = .TRUE.
      91         472 :          needs%norm_drho_spin = .TRUE.
      92             :       END IF
      93         473 :       IF (PRESENT(max_deriv)) max_deriv = 1
      94         473 :    END SUBROUTINE optx_lsd_info
      95             : 
      96             : ! **************************************************************************************************
      97             : !> \brief evaluates the optx functional for lda
      98             : !> \param rho_set the density where you want to evaluate the functional
      99             : !> \param deriv_set place where to store the functional derivatives (they are
     100             : !>        added to the derivatives)
     101             : !> \param grad_deriv degree of the derivative that should be evaluated,
     102             : !>        if positive all the derivatives up to the given degree are evaluated,
     103             : !>        if negative only the given degree is calculated
     104             : !> \param optx_params input parameter (scaling)
     105             : !> \par History
     106             : !>      01.2007 added scaling [Manuel Guidon]
     107             : !> \author Joost
     108             : ! **************************************************************************************************
     109        2820 :    SUBROUTINE optx_lda_eval(rho_set, deriv_set, grad_deriv, optx_params)
     110             :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     111             :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     112             :       INTEGER, INTENT(in)                                :: grad_deriv
     113             :       TYPE(section_vals_type), POINTER                   :: optx_params
     114             : 
     115             :       INTEGER                                            :: npoints
     116             :       INTEGER, DIMENSION(2, 3)                           :: bo
     117             :       REAL(kind=dp)                                      :: a1, a2, epsilon_drho, epsilon_rho, gam, &
     118             :                                                             sx
     119             :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     120         564 :          POINTER                                         :: e_0, e_ndrho, e_rho, norm_drho, rho
     121             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     122             : 
     123         564 :       NULLIFY (e_0, e_ndrho, e_rho, norm_drho, rho)
     124             : 
     125         564 :       CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
     126         564 :       CALL section_vals_val_get(optx_params, "a1", r_val=a1)
     127         564 :       CALL section_vals_val_get(optx_params, "a2", r_val=a2)
     128         564 :       CALL section_vals_val_get(optx_params, "gamma", r_val=gam)
     129             : 
     130             :       CALL xc_rho_set_get(rho_set, rho=rho, &
     131             :                           norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
     132         564 :                           drho_cutoff=epsilon_drho)
     133         564 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     134             : 
     135             :       deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     136         564 :                                       allocate_deriv=.TRUE.)
     137         564 :       CALL xc_derivative_get(deriv, deriv_data=e_0)
     138             :       deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     139         564 :                                       allocate_deriv=.TRUE.)
     140         564 :       CALL xc_derivative_get(deriv, deriv_data=e_rho)
     141             :       deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     142         564 :                                       allocate_deriv=.TRUE.)
     143         564 :       CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     144         564 :       IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
     145           0 :          CPABORT("derivatives bigger than 1 not implemented")
     146             :       END IF
     147             : 
     148             :       CALL optx_lda_calc(rho=rho, norm_drho=norm_drho, &
     149             :                          e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
     150             :                          npoints=npoints, epsilon_rho=epsilon_rho, &
     151             :                          epsilon_drho=epsilon_drho, sx=sx, &
     152         564 :                          a1=a1, a2=a2, gam=gam)
     153         564 :    END SUBROUTINE optx_lda_eval
     154             : 
     155             : ! **************************************************************************************************
     156             : !> \brief evaluates the optx functional for lsd
     157             : !> \param rho_set the density where you want to evaluate the functional
     158             : !> \param deriv_set place where to store the functional derivatives (they are
     159             : !>        added to the derivatives)
     160             : !> \param grad_deriv degree of the derivative that should be evaluated,
     161             : !>        if positive all the derivatives up to the given degree are evaluated,
     162             : !>        if negative only the given degree is calculated
     163             : !> \param optx_params input parameter (scaling)
     164             : !> \par History
     165             : !>      01.2007 added scaling [Manuel Guidon]
     166             : !> \author Joost
     167             : ! **************************************************************************************************
     168        2808 :    SUBROUTINE optx_lsd_eval(rho_set, deriv_set, grad_deriv, optx_params)
     169             :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     170             :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     171             :       INTEGER, INTENT(in)                                :: grad_deriv
     172             :       TYPE(section_vals_type), POINTER                   :: optx_params
     173             : 
     174             :       INTEGER                                            :: ispin, npoints
     175             :       INTEGER, DIMENSION(2, 3)                           :: bo
     176             :       REAL(kind=dp)                                      :: a1, a2, epsilon_drho, epsilon_rho, gam, &
     177             :                                                             sx
     178             :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     179         468 :          POINTER                                         :: e_0
     180        5616 :       TYPE(cp_3d_r_cp_type), DIMENSION(2)                :: e_ndrho, e_rho, ndrho, rho
     181             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     182             : 
     183         468 :       NULLIFY (e_0)
     184        1404 :       DO ispin = 1, 2
     185         936 :          NULLIFY (e_rho(ispin)%array)
     186         936 :          NULLIFY (e_ndrho(ispin)%array)
     187         936 :          NULLIFY (rho(ispin)%array)
     188        1404 :          NULLIFY (ndrho(ispin)%array)
     189             :       END DO
     190             : 
     191         468 :       CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
     192         468 :       CALL section_vals_val_get(optx_params, "a1", r_val=a1)
     193         468 :       CALL section_vals_val_get(optx_params, "a2", r_val=a2)
     194         468 :       CALL section_vals_val_get(optx_params, "gamma", r_val=gam)
     195             : 
     196             :       CALL xc_rho_set_get(rho_set, rhoa=rho(1)%array, rhob=rho(2)%array, &
     197             :                           norm_drhoa=ndrho(1)%array, &
     198             :                           norm_drhob=ndrho(2)%array, rho_cutoff=epsilon_rho, &
     199         468 :                           drho_cutoff=epsilon_drho, local_bounds=bo)
     200         468 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     201             : 
     202             :       deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     203         468 :                                       allocate_deriv=.TRUE.)
     204         468 :       CALL xc_derivative_get(deriv, deriv_data=e_0)
     205             :       deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
     206         468 :                                       allocate_deriv=.TRUE.)
     207         468 :       CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
     208             :       deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
     209         468 :                                       allocate_deriv=.TRUE.)
     210         468 :       CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
     211             : 
     212             :       deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
     213         468 :                                       allocate_deriv=.TRUE.)
     214         468 :       CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
     215             :       deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
     216         468 :                                       allocate_deriv=.TRUE.)
     217         468 :       CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
     218             : 
     219         468 :       IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
     220           0 :          CPABORT("derivatives bigger than 1 not implemented")
     221             :       END IF
     222        1404 :       DO ispin = 1, 2
     223             :          CALL optx_lsd_calc(rho=rho(ispin)%array, norm_drho=ndrho(ispin)%array, &
     224             :                             e_0=e_0, e_rho=e_rho(ispin)%array, e_ndrho=e_ndrho(ispin)%array, &
     225             :                             npoints=npoints, epsilon_rho=epsilon_rho, &
     226             :                             epsilon_drho=epsilon_drho, sx=sx, &
     227        1404 :                             a1=a1, a2=a2, gam=gam)
     228             :       END DO
     229         468 :    END SUBROUTINE optx_lsd_eval
     230             : 
     231             : ! **************************************************************************************************
     232             : !> \brief optx exchange functional
     233             : !> \param rho the full density
     234             : !> \param norm_drho the norm of the gradient of the full density
     235             : !> \param e_0 the value of the functional in that point
     236             : !> \param e_rho the derivative of the functional wrt. rho
     237             : !> \param e_ndrho the derivative of the functional wrt. norm_drho
     238             : !> \param epsilon_rho the cutoff on rho
     239             : !> \param epsilon_drho ...
     240             : !> \param npoints ...
     241             : !> \param sx scaling-parameter for exchange
     242             : !> \param a1 a1 coefficient of the OPTX functional
     243             : !> \param a2 a2 coefficient of the OPTX functional
     244             : !> \param gam gamma coefficient of the OPTX functional
     245             : !> \par History
     246             : !>      01.2007 added scaling [Manuel Guidon]
     247             : !> \author Joost VandeVondele
     248             : ! **************************************************************************************************
     249         564 :    SUBROUTINE optx_lda_calc(rho, norm_drho, e_0, e_rho, e_ndrho, &
     250             :                             epsilon_rho, epsilon_drho, npoints, sx, &
     251             :                             a1, a2, gam)
     252             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, norm_drho
     253             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_0, e_rho, e_ndrho
     254             :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho, epsilon_drho
     255             :       INTEGER, INTENT(in)                                :: npoints
     256             :       REAL(kind=dp), INTENT(in)                          :: sx, a1, a2, gam
     257             : 
     258             :       REAL(KIND=dp), PARAMETER                           :: cx = 0.930525736349100_dp, &
     259             :                                                             o43 = 4.0_dp/3.0_dp
     260             : 
     261             :       INTEGER                                            :: ii
     262             :       REAL(KIND=dp)                                      :: denom, ex, gamxsxs, myndrho, myrho, &
     263             :                                                             rho43, tmp, xs
     264             : 
     265             : !$OMP     PARALLEL DO DEFAULT (NONE) &
     266             : !$OMP                 SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
     267             : !$OMP                 SHARED(epsilon_rho, epsilon_drho, sx, npoints) &
     268             : !$OMP                 SHARED(a1, a2, gam) &
     269             : !$OMP                 PRIVATE(ii, myrho, myndrho, rho43, xs, gamxsxs) &
     270         564 : !$OMP                 PRIVATE(denom, ex, tmp)
     271             : 
     272             :       DO ii = 1, npoints
     273             :          ! we get the full density and need spin parts -> 0.5
     274             :          myrho = 0.5_dp*rho(ii)
     275             :          myndrho = 0.5_dp*MAX(norm_drho(ii), epsilon_drho)
     276             :          IF (myrho > 0.5_dp*epsilon_rho) THEN
     277             :             rho43 = myrho**o43
     278             :             xs = (myndrho/rho43)
     279             :             gamxsxs = gam*xs*xs
     280             :             denom = 1.0_dp/(1.0_dp + gamxsxs)
     281             :             ex = rho43*(a1*cx + a2*(gamxsxs*denom)**2)
     282             :             ! 2.0 for both spins
     283             :             e_0(ii) = e_0(ii) - (2.0_dp*ex)*sx
     284             :             tmp = rho43*2.0_dp*a2*gamxsxs*denom**2*(1.0_dp - gamxsxs*denom)
     285             :             ! derive e_0 wrt to rho (full) and ndrho (also full)
     286             :             e_rho(ii) = e_rho(ii) - ((o43*ex + tmp*gamxsxs*(-2.0_dp*o43))/myrho)*sx
     287             :             e_ndrho(ii) = e_ndrho(ii) - ((tmp*gam*2.0_dp*myndrho/rho43**2))*sx
     288             :          END IF
     289             :       END DO
     290             : 
     291             : !$OMP     END PARALLEL DO
     292             : 
     293         564 :    END SUBROUTINE optx_lda_calc
     294             : 
     295             : ! **************************************************************************************************
     296             : !> \brief optx exchange functional
     297             : !> \param rho the *spin* density
     298             : !> \param norm_drho the norm of the gradient of the *spin* density
     299             : !> \param e_0 the value of the functional in that point
     300             : !> \param e_rho the derivative of the functional wrt. rho
     301             : !> \param e_ndrho the derivative of the functional wrt. norm_drho
     302             : !> \param epsilon_rho the cutoff on rho
     303             : !> \param epsilon_drho ...
     304             : !> \param npoints ...
     305             : !> \param sx scaling parameter for exchange
     306             : !> \param a1 a1 coefficient of the OPTX functional
     307             : !> \param a2 a2 coefficient of the OPTX functional
     308             : !> \param gam gamma coefficient of the OPTX functional
     309             : !> \par History
     310             : !>      01.2007 added scaling [Manuel Guidon]
     311             : !> \author Joost VandeVondele
     312             : ! **************************************************************************************************
     313         936 :    SUBROUTINE optx_lsd_calc(rho, norm_drho, e_0, e_rho, e_ndrho, &
     314             :                             epsilon_rho, epsilon_drho, npoints, sx, &
     315             :                             a1, a2, gam)
     316             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, norm_drho
     317             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_0, e_rho, e_ndrho
     318             :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho, epsilon_drho
     319             :       INTEGER, INTENT(in)                                :: npoints
     320             :       REAL(kind=dp), INTENT(in)                          :: sx, a1, a2, gam
     321             : 
     322             :       REAL(KIND=dp), PARAMETER                           :: cx = 0.930525736349100_dp, &
     323             :                                                             o43 = 4.0_dp/3.0_dp
     324             : 
     325             :       INTEGER                                            :: ii
     326             :       REAL(KIND=dp)                                      :: denom, ex, gamxsxs, myndrho, myrho, &
     327             :                                                             rho43, tmp, xs
     328             : 
     329             : !$OMP     PARALLEL DO DEFAULT(NONE) &
     330             : !$OMP                 SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
     331             : !$OMP                 SHARED(epsilon_rho, epsilon_drho, npoints, sx) &
     332             : !$OMP                 SHARED(a1, a2, gam) &
     333             : !$OMP                 PRIVATE(ii, denom, ex, gamxsxs, myndrho, myrho) &
     334         936 : !$OMP                 PRIVATE(rho43, tmp, xs)
     335             : 
     336             :       DO ii = 1, npoints
     337             :          ! we do have the spin density already
     338             :          myrho = rho(ii)
     339             :          myndrho = MAX(norm_drho(ii), epsilon_drho)
     340             :          IF (myrho > epsilon_rho) THEN
     341             :             rho43 = myrho**o43
     342             :             xs = (myndrho/rho43)
     343             :             gamxsxs = gam*xs*xs
     344             :             denom = 1.0_dp/(1.0_dp + gamxsxs)
     345             :             ex = rho43*(a1*cx + a2*(gamxsxs*denom)**2)
     346             :             ! for a single spin
     347             :             e_0(ii) = e_0(ii) - ex*sx
     348             :             tmp = rho43*2.0_dp*a2*gamxsxs*denom**2*(1.0_dp - gamxsxs*denom)
     349             :             ! derive e_0 wrt to rho and ndrho
     350             :             e_rho(ii) = e_rho(ii) - ((o43*ex + tmp*gamxsxs*(-2.0_dp*o43))/myrho)*sx
     351             :             e_ndrho(ii) = e_ndrho(ii) - ((tmp*gam*2.0_dp*myndrho/rho43**2))*sx
     352             :          END IF
     353             :       END DO
     354             : 
     355             : !$OMP     END PARALLEL DO
     356             : 
     357         936 :    END SUBROUTINE optx_lsd_calc
     358             : 
     359             : END MODULE xc_optx