LCOV - code coverage report
Current view: top level - src - qs_fxc.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 56.3 % 238 134
Test Date: 2025-07-25 12:55:17 Functions: 80.0 % 5 4

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief  https://en.wikipedia.org/wiki/Finite_difference_coefficient
      10              : !---------------------------------------------------------------------------------------------------
      11              : !Derivative    Accuracy            4       3       2       1       0       1       2      3       4
      12              : !---------------------------------------------------------------------------------------------------
      13              : !    1             2                                    -1/2       0     1/2
      14              : !                  4                            1/12    -2/3       0     2/3   -1/12
      15              : !                  6                   -1/60    3/20    -3/4       0     3/4   -3/20   1/60
      16              : !                  8           1/280  -4/105     1/5    -4/5       0     4/5    -1/5  4/105  -1/280
      17              : !---------------------------------------------------------------------------------------------------
      18              : !    2             2                                       1      -2       1
      19              : !                  4                           -1/12     4/3    -5/2     4/3   -1/12
      20              : !                  6                    1/90   -3/20     3/2  -49/18     3/2   -3/20   1/90
      21              : !                  8          -1/560   8/315    -1/5     8/5 -205/72     8/5    -1/5  8/315  -1/560
      22              : !---------------------------------------------------------------------------------------------------
      23              : !> \par History
      24              : !>     init 17.03.2020
      25              : !> \author JGH
      26              : ! **************************************************************************************************
      27              : MODULE qs_fxc
      28              : 
      29              :    USE cp_control_types,                ONLY: dft_control_type
      30              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      31              :                                               section_vals_type
      32              :    USE kinds,                           ONLY: dp
      33              :    USE pw_env_types,                    ONLY: pw_env_get,&
      34              :                                               pw_env_type
      35              :    USE pw_methods,                      ONLY: pw_axpy,&
      36              :                                               pw_scale,&
      37              :                                               pw_zero
      38              :    USE pw_pool_types,                   ONLY: pw_pool_type
      39              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      40              :                                               pw_r3d_rs_type
      41              :    USE qs_ks_types,                     ONLY: get_ks_env,&
      42              :                                               qs_ks_env_type
      43              :    USE qs_rho_methods,                  ONLY: qs_rho_copy,&
      44              :                                               qs_rho_scale_and_add
      45              :    USE qs_rho_types,                    ONLY: qs_rho_create,&
      46              :                                               qs_rho_get,&
      47              :                                               qs_rho_release,&
      48              :                                               qs_rho_type
      49              :    USE qs_vxc,                          ONLY: qs_vxc_create
      50              :    USE xc,                              ONLY: xc_calc_2nd_deriv,&
      51              :                                               xc_prep_2nd_deriv
      52              :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      53              :                                               xc_dset_release
      54              :    USE xc_derivatives,                  ONLY: xc_functionals_get_needs
      55              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      56              :    USE xc_rho_set_types,                ONLY: xc_rho_set_release,&
      57              :                                               xc_rho_set_type
      58              : #include "./base/base_uses.f90"
      59              : 
      60              :    IMPLICIT NONE
      61              : 
      62              :    PRIVATE
      63              : 
      64              :    ! *** Public subroutines ***
      65              :    PUBLIC :: qs_fxc_fdiff, qs_fxc_analytic, qs_fgxc_gdiff, qs_fgxc_create, qs_fgxc_release
      66              : 
      67              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fxc'
      68              : 
      69              : ! **************************************************************************************************
      70              : 
      71              : CONTAINS
      72              : 
      73              : ! **************************************************************************************************
      74              : !> \brief ...
      75              : !> \param rho0 ...
      76              : !> \param rho1_r ...
      77              : !> \param tau1_r ...
      78              : !> \param xc_section ...
      79              : !> \param auxbas_pw_pool ...
      80              : !> \param is_triplet ...
      81              : !> \param v_xc ...
      82              : !> \param v_xc_tau ...
      83              : ! **************************************************************************************************
      84        16380 :    SUBROUTINE qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau)
      85              : 
      86              :       TYPE(qs_rho_type), POINTER                         :: rho0
      87              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho1_r, tau1_r
      88              :       TYPE(section_vals_type), POINTER                   :: xc_section
      89              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
      90              :       LOGICAL, INTENT(IN)                                :: is_triplet
      91              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_xc, v_xc_tau
      92              : 
      93              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_fxc_analytic'
      94              : 
      95              :       INTEGER                                            :: handle, nspins
      96              :       INTEGER, DIMENSION(2, 3)                           :: bo
      97              :       LOGICAL                                            :: lsd
      98              :       REAL(KIND=dp)                                      :: fac
      99        16380 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho0_g, rho1_g
     100        16380 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho0_r, tau0_r
     101              :       TYPE(section_vals_type), POINTER                   :: xc_fun_section
     102              :       TYPE(xc_derivative_set_type)                       :: deriv_set
     103              :       TYPE(xc_rho_cflags_type)                           :: needs
     104              :       TYPE(xc_rho_set_type)                              :: rho0_set
     105              : 
     106         8190 :       CALL timeset(routineN, handle)
     107              : 
     108         8190 :       CPASSERT(.NOT. ASSOCIATED(v_xc))
     109         8190 :       CPASSERT(.NOT. ASSOCIATED(v_xc_tau))
     110              : 
     111         8190 :       CALL qs_rho_get(rho0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
     112         8190 :       nspins = SIZE(rho0_r)
     113              : 
     114         8190 :       lsd = (nspins == 2)
     115              :       fac = 0._dp
     116         8190 :       IF (is_triplet .AND. nspins == 1) fac = -1.0_dp
     117              : 
     118         8190 :       NULLIFY (rho1_g)
     119        81900 :       bo = rho1_r(1)%pw_grid%bounds_local
     120         8190 :       xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
     121         8190 :       needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
     122              :       ! calculate the arguments needed by the functionals
     123         8190 :       CALL xc_prep_2nd_deriv(deriv_set, rho0_set, rho0_r, auxbas_pw_pool, xc_section=xc_section, tau_r=tau0_r)
     124              :       CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho0_set, rho1_r, rho1_g, tau1_r, &
     125         8190 :                              auxbas_pw_pool, xc_section=xc_section, gapw=.FALSE., do_triplet=is_triplet)
     126         8190 :       CALL xc_dset_release(deriv_set)
     127         8190 :       CALL xc_rho_set_release(rho0_set)
     128              : 
     129         8190 :       CALL timestop(handle)
     130              : 
     131       180180 :    END SUBROUTINE qs_fxc_analytic
     132              : 
     133              : ! **************************************************************************************************
     134              : !> \brief ...
     135              : !> \param ks_env ...
     136              : !> \param rho0_struct ...
     137              : !> \param rho1_struct ...
     138              : !> \param xc_section ...
     139              : !> \param accuracy ...
     140              : !> \param is_triplet ...
     141              : !> \param fxc_rho ...
     142              : !> \param fxc_tau ...
     143              : ! **************************************************************************************************
     144         1908 :    SUBROUTINE qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
     145              :                            fxc_rho, fxc_tau)
     146              : 
     147              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     148              :       TYPE(qs_rho_type), POINTER                         :: rho0_struct, rho1_struct
     149              :       TYPE(section_vals_type), POINTER                   :: xc_section
     150              :       INTEGER, INTENT(IN)                                :: accuracy
     151              :       LOGICAL, INTENT(IN)                                :: is_triplet
     152              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: fxc_rho, fxc_tau
     153              : 
     154              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_fxc_fdiff'
     155              :       REAL(KIND=dp), PARAMETER                           :: epsrho = 5.e-4_dp
     156              : 
     157              :       INTEGER                                            :: handle, ispin, istep, nspins, nstep
     158              :       REAL(KIND=dp)                                      :: alpha, beta, exc, oeps1
     159              :       REAL(KIND=dp), DIMENSION(-4:4)                     :: ak
     160              :       TYPE(dft_control_type), POINTER                    :: dft_control
     161              :       TYPE(pw_env_type), POINTER                         :: pw_env
     162              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     163         1908 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_tau_rspace, vxc00
     164              :       TYPE(qs_rho_type), POINTER                         :: rhoin
     165              : 
     166         1908 :       CALL timeset(routineN, handle)
     167              : 
     168         1908 :       CPASSERT(.NOT. ASSOCIATED(fxc_rho))
     169         1908 :       CPASSERT(.NOT. ASSOCIATED(fxc_tau))
     170         1908 :       CPASSERT(ASSOCIATED(rho0_struct))
     171         1908 :       CPASSERT(ASSOCIATED(rho1_struct))
     172              : 
     173         1908 :       ak = 0.0_dp
     174         1908 :       SELECT CASE (accuracy)
     175              :       CASE (:4)
     176            0 :          nstep = 2
     177            0 :          ak(-2:2) = (/1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp/)/12.0_dp
     178              :       CASE (5:7)
     179        15264 :          nstep = 3
     180        15264 :          ak(-3:3) = (/-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp/)/60.0_dp
     181              :       CASE (8:)
     182            0 :          nstep = 4
     183              :          ak(-4:4) = (/1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
     184         1908 :                       224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp/)/280.0_dp
     185              :       END SELECT
     186              : 
     187         1908 :       CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
     188         1908 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     189              : 
     190         1908 :       nspins = dft_control%nspins
     191         1908 :       exc = 0.0_dp
     192              : 
     193        15264 :       DO istep = -nstep, nstep
     194              : 
     195        15264 :          IF (ak(istep) /= 0.0_dp) THEN
     196        11448 :             alpha = 1.0_dp
     197        11448 :             beta = REAL(istep, KIND=dp)*epsrho
     198              :             NULLIFY (rhoin)
     199        11448 :             ALLOCATE (rhoin)
     200        11448 :             CALL qs_rho_create(rhoin)
     201        11448 :             NULLIFY (vxc00, v_tau_rspace)
     202        11448 :             IF (is_triplet) THEN
     203          924 :                CPASSERT(nspins == 1)
     204              :                ! rhoin = (0.5 rho0, 0.5 rho0)
     205          924 :                CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
     206              :                ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
     207          924 :                CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
     208              :                CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
     209          924 :                                   vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
     210          924 :                CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
     211          924 :                IF (ASSOCIATED(v_tau_rspace)) CALL pw_axpy(v_tau_rspace(2), v_tau_rspace(1), -1.0_dp)
     212              :             ELSE
     213        10524 :                CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
     214        10524 :                CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
     215              :                CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
     216        10524 :                                   vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
     217              :             END IF
     218        11448 :             CALL qs_rho_release(rhoin)
     219        11448 :             DEALLOCATE (rhoin)
     220        11448 :             IF (.NOT. ASSOCIATED(fxc_rho)) THEN
     221         7932 :                ALLOCATE (fxc_rho(nspins))
     222         4116 :                DO ispin = 1, nspins
     223         2208 :                   CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
     224         4116 :                   CALL pw_zero(fxc_rho(ispin))
     225              :                END DO
     226              :             END IF
     227        24696 :             DO ispin = 1, nspins
     228        24696 :                CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
     229              :             END DO
     230        25620 :             DO ispin = 1, SIZE(vxc00)
     231        25620 :                CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
     232              :             END DO
     233        11448 :             DEALLOCATE (vxc00)
     234        11448 :             IF (ASSOCIATED(v_tau_rspace)) THEN
     235            0 :                IF (.NOT. ASSOCIATED(fxc_tau)) THEN
     236            0 :                   ALLOCATE (fxc_tau(nspins))
     237            0 :                   DO ispin = 1, nspins
     238            0 :                      CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
     239            0 :                      CALL pw_zero(fxc_tau(ispin))
     240              :                   END DO
     241              :                END IF
     242            0 :                DO ispin = 1, nspins
     243            0 :                   CALL pw_axpy(v_tau_rspace(ispin), fxc_tau(ispin), ak(istep))
     244              :                END DO
     245            0 :                DO ispin = 1, SIZE(v_tau_rspace)
     246            0 :                   CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
     247              :                END DO
     248            0 :                DEALLOCATE (v_tau_rspace)
     249              :             END IF
     250              :          END IF
     251              : 
     252              :       END DO
     253              : 
     254         1908 :       oeps1 = 1.0_dp/epsrho
     255         4116 :       DO ispin = 1, nspins
     256         4116 :          CALL pw_scale(fxc_rho(ispin), oeps1)
     257              :       END DO
     258         1908 :       IF (ASSOCIATED(fxc_tau)) THEN
     259            0 :          DO ispin = 1, nspins
     260            0 :             CALL pw_scale(fxc_tau(ispin), oeps1)
     261              :          END DO
     262              :       END IF
     263              : 
     264         1908 :       CALL timestop(handle)
     265              : 
     266         1908 :    END SUBROUTINE qs_fxc_fdiff
     267              : 
     268              : ! **************************************************************************************************
     269              : !> \brief ...
     270              : !> \param ks_env ...
     271              : !> \param rho0_struct ...
     272              : !> \param rho1_struct ...
     273              : !> \param xc_section ...
     274              : !> \param accuracy ...
     275              : !> \param epsrho ...
     276              : !> \param is_triplet ...
     277              : !> \param fxc_rho ...
     278              : !> \param fxc_tau ...
     279              : !> \param gxc_rho ...
     280              : !> \param gxc_tau ...
     281              : ! **************************************************************************************************
     282          264 :    SUBROUTINE qs_fgxc_gdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, &
     283              :                             is_triplet, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     284              : 
     285              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     286              :       TYPE(qs_rho_type), POINTER                         :: rho0_struct, rho1_struct
     287              :       TYPE(section_vals_type), POINTER                   :: xc_section
     288              :       INTEGER, INTENT(IN)                                :: accuracy
     289              :       REAL(KIND=dp), INTENT(IN)                          :: epsrho
     290              :       LOGICAL, INTENT(IN)                                :: is_triplet
     291              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
     292              : 
     293              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_fgxc_gdiff'
     294              : 
     295              :       INTEGER                                            :: handle, ispin, istep, nspins, nstep
     296              :       REAL(KIND=dp)                                      :: alpha, beta, exc, oeps1
     297              :       REAL(KIND=dp), DIMENSION(-4:4)                     :: ak
     298              :       TYPE(dft_control_type), POINTER                    :: dft_control
     299              :       TYPE(pw_env_type), POINTER                         :: pw_env
     300              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     301          264 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_tau_rspace, vxc00
     302              :       TYPE(qs_rho_type), POINTER                         :: rhoin
     303              : 
     304          264 :       CALL timeset(routineN, handle)
     305              : 
     306          264 :       CPASSERT(.NOT. ASSOCIATED(fxc_rho))
     307          264 :       CPASSERT(.NOT. ASSOCIATED(fxc_tau))
     308          264 :       CPASSERT(.NOT. ASSOCIATED(gxc_rho))
     309          264 :       CPASSERT(.NOT. ASSOCIATED(gxc_tau))
     310          264 :       CPASSERT(ASSOCIATED(rho0_struct))
     311          264 :       CPASSERT(ASSOCIATED(rho1_struct))
     312              : 
     313          264 :       ak = 0.0_dp
     314          264 :       SELECT CASE (accuracy)
     315              :       CASE (:4)
     316            0 :          nstep = 2
     317            0 :          ak(-2:2) = (/1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp/)/12.0_dp
     318              :       CASE (5:7)
     319         2112 :          nstep = 3
     320         2112 :          ak(-3:3) = (/-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp/)/60.0_dp
     321              :       CASE (8:)
     322            0 :          nstep = 4
     323              :          ak(-4:4) = (/1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
     324          264 :                       224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp/)/280.0_dp
     325              :       END SELECT
     326              : 
     327          264 :       CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
     328          264 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     329              : 
     330          264 :       nspins = dft_control%nspins
     331          264 :       exc = 0.0_dp
     332              : 
     333              :       CALL qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
     334          264 :                         fxc_rho, fxc_tau)
     335              : 
     336         2112 :       DO istep = -nstep, nstep
     337              : 
     338         2112 :          IF (ak(istep) /= 0.0_dp) THEN
     339         1584 :             alpha = 1.0_dp
     340         1584 :             beta = REAL(istep, KIND=dp)*epsrho
     341              :             NULLIFY (rhoin)
     342         1584 :             ALLOCATE (rhoin)
     343         1584 :             CALL qs_rho_create(rhoin)
     344         1584 :             NULLIFY (vxc00, v_tau_rspace)
     345         1584 :             CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
     346         1584 :             CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
     347              :             CALL qs_fxc_fdiff(ks_env=ks_env, rho0_struct=rhoin, rho1_struct=rho1_struct, &
     348              :                               xc_section=xc_section, accuracy=accuracy, is_triplet=is_triplet, &
     349         1584 :                               fxc_rho=vxc00, fxc_tau=v_tau_rspace)
     350         1584 :             CALL qs_rho_release(rhoin)
     351         1584 :             DEALLOCATE (rhoin)
     352         1584 :             IF (.NOT. ASSOCIATED(gxc_rho)) THEN
     353         1092 :                ALLOCATE (gxc_rho(nspins))
     354          564 :                DO ispin = 1, nspins
     355          300 :                   CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
     356          564 :                   CALL pw_zero(gxc_rho(ispin))
     357              :                END DO
     358              :             END IF
     359         3384 :             DO ispin = 1, nspins
     360         3384 :                CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
     361              :             END DO
     362         3384 :             DO ispin = 1, SIZE(vxc00)
     363         3384 :                CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
     364              :             END DO
     365         1584 :             DEALLOCATE (vxc00)
     366         1584 :             IF (ASSOCIATED(v_tau_rspace)) THEN
     367            0 :                IF (.NOT. ASSOCIATED(gxc_tau)) THEN
     368            0 :                   ALLOCATE (gxc_tau(nspins))
     369            0 :                   DO ispin = 1, nspins
     370            0 :                      CALL auxbas_pw_pool%create_pw(gxc_tau(ispin))
     371            0 :                      CALL pw_zero(gxc_tau(ispin))
     372              :                   END DO
     373              :                END IF
     374            0 :                DO ispin = 1, nspins
     375            0 :                   CALL pw_axpy(v_tau_rspace(ispin), gxc_tau(ispin), ak(istep))
     376              :                END DO
     377            0 :                DO ispin = 1, SIZE(v_tau_rspace)
     378            0 :                   CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
     379              :                END DO
     380            0 :                DEALLOCATE (v_tau_rspace)
     381              :             END IF
     382              :          END IF
     383              : 
     384              :       END DO
     385              : 
     386          264 :       oeps1 = 1.0_dp/epsrho
     387          564 :       DO ispin = 1, nspins
     388          564 :          CALL pw_scale(gxc_rho(ispin), oeps1)
     389              :       END DO
     390          264 :       IF (ASSOCIATED(gxc_tau)) THEN
     391            0 :          DO ispin = 1, nspins
     392            0 :             CALL pw_scale(gxc_tau(ispin), oeps1)
     393              :          END DO
     394              :       END IF
     395              : 
     396          264 :       CALL timestop(handle)
     397              : 
     398          264 :    END SUBROUTINE qs_fgxc_gdiff
     399              : 
     400              : ! **************************************************************************************************
     401              : !> \brief ...
     402              : !> \param ks_env ...
     403              : !> \param rho0_struct ...
     404              : !> \param rho1_struct ...
     405              : !> \param xc_section ...
     406              : !> \param accuracy ...
     407              : !> \param is_triplet ...
     408              : !> \param fxc_rho ...
     409              : !> \param fxc_tau ...
     410              : !> \param gxc_rho ...
     411              : !> \param gxc_tau ...
     412              : ! **************************************************************************************************
     413            0 :    SUBROUTINE qs_fgxc_create(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
     414              :                              fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     415              : 
     416              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     417              :       TYPE(qs_rho_type), POINTER                         :: rho0_struct, rho1_struct
     418              :       TYPE(section_vals_type), POINTER                   :: xc_section
     419              :       INTEGER, INTENT(IN)                                :: accuracy
     420              :       LOGICAL, INTENT(IN)                                :: is_triplet
     421              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
     422              : 
     423              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_fgxc_create'
     424              :       REAL(KIND=dp), PARAMETER                           :: epsrho = 5.e-4_dp
     425              : 
     426              :       INTEGER                                            :: handle, ispin, istep, nspins, nstep
     427              :       REAL(KIND=dp)                                      :: alpha, beta, exc, oeps1, oeps2
     428              :       REAL(KIND=dp), DIMENSION(-4:4)                     :: ak, bl
     429              :       TYPE(dft_control_type), POINTER                    :: dft_control
     430              :       TYPE(pw_env_type), POINTER                         :: pw_env
     431              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     432            0 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_tau_rspace, vxc00
     433              :       TYPE(qs_rho_type), POINTER                         :: rhoin
     434              : 
     435            0 :       CALL timeset(routineN, handle)
     436              : 
     437            0 :       CPASSERT(.NOT. ASSOCIATED(fxc_rho))
     438            0 :       CPASSERT(.NOT. ASSOCIATED(fxc_tau))
     439            0 :       CPASSERT(.NOT. ASSOCIATED(gxc_rho))
     440            0 :       CPASSERT(.NOT. ASSOCIATED(gxc_tau))
     441            0 :       CPASSERT(ASSOCIATED(rho0_struct))
     442            0 :       CPASSERT(ASSOCIATED(rho1_struct))
     443              : 
     444            0 :       ak = 0.0_dp
     445            0 :       bl = 0.0_dp
     446            0 :       SELECT CASE (accuracy)
     447              :       CASE (:4)
     448            0 :          nstep = 2
     449            0 :          ak(-2:2) = (/1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp/)/12.0_dp
     450            0 :          bl(-2:2) = (/-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp/)/12.0_dp
     451              :       CASE (5:7)
     452            0 :          nstep = 3
     453            0 :          ak(-3:3) = (/-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp/)/60.0_dp
     454            0 :          bl(-3:3) = (/2.0_dp, -27.0_dp, 270.0_dp, -490.0_dp, 270.0_dp, -27.0_dp, 2.0_dp/)/180.0_dp
     455              :       CASE (8:)
     456            0 :          nstep = 4
     457              :          ak(-4:4) = (/1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
     458            0 :                       224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp/)/280.0_dp
     459              :          bl(-4:4) = (/-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
     460            0 :                       896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp/)/560.0_dp
     461              :       END SELECT
     462              : 
     463            0 :       CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
     464            0 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     465              : 
     466            0 :       nspins = dft_control%nspins
     467            0 :       exc = 0.0_dp
     468              : 
     469            0 :       DO istep = -nstep, nstep
     470              : 
     471            0 :          alpha = 1.0_dp
     472            0 :          beta = REAL(istep, KIND=dp)*epsrho
     473              :          NULLIFY (rhoin)
     474            0 :          ALLOCATE (rhoin)
     475            0 :          CALL qs_rho_create(rhoin)
     476            0 :          NULLIFY (vxc00, v_tau_rspace)
     477            0 :          IF (is_triplet) THEN
     478            0 :             CPASSERT(nspins == 1)
     479              :             ! rhoin = (0.5 rho0, 0.5 rho0)
     480            0 :             CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
     481              :             ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
     482            0 :             CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
     483              :             CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
     484            0 :                                vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
     485            0 :             CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
     486              :          ELSE
     487            0 :             CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
     488            0 :             CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
     489              :             CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
     490            0 :                                vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
     491              :          END IF
     492            0 :          CALL qs_rho_release(rhoin)
     493            0 :          DEALLOCATE (rhoin)
     494            0 :          IF (.NOT. ASSOCIATED(fxc_rho)) THEN
     495            0 :             ALLOCATE (fxc_rho(nspins))
     496            0 :             DO ispin = 1, nspins
     497            0 :                CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
     498            0 :                CALL pw_zero(fxc_rho(ispin))
     499              :             END DO
     500              :          END IF
     501            0 :          IF (.NOT. ASSOCIATED(gxc_rho)) THEN
     502            0 :             ALLOCATE (gxc_rho(nspins))
     503            0 :             DO ispin = 1, nspins
     504            0 :                CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
     505            0 :                CALL pw_zero(gxc_rho(ispin))
     506              :             END DO
     507              :          END IF
     508            0 :          CPASSERT(.NOT. ASSOCIATED(v_tau_rspace))
     509            0 :          DO ispin = 1, nspins
     510            0 :             IF (ak(istep) /= 0.0_dp) THEN
     511            0 :                CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
     512              :             END IF
     513            0 :             IF (bl(istep) /= 0.0_dp) THEN
     514            0 :                CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), bl(istep))
     515              :             END IF
     516              :          END DO
     517            0 :          DO ispin = 1, SIZE(vxc00)
     518            0 :             CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
     519              :          END DO
     520            0 :          DEALLOCATE (vxc00)
     521              : 
     522              :       END DO
     523              : 
     524            0 :       oeps1 = 1.0_dp/epsrho
     525            0 :       oeps2 = 1.0_dp/(epsrho**2)
     526            0 :       DO ispin = 1, nspins
     527            0 :          CALL pw_scale(fxc_rho(ispin), oeps1)
     528            0 :          CALL pw_scale(gxc_rho(ispin), oeps2)
     529              :       END DO
     530              : 
     531            0 :       CALL timestop(handle)
     532              : 
     533            0 :    END SUBROUTINE qs_fgxc_create
     534              : 
     535              : ! **************************************************************************************************
     536              : !> \brief ...
     537              : !> \param ks_env ...
     538              : !> \param fxc_rho ...
     539              : !> \param fxc_tau ...
     540              : !> \param gxc_rho ...
     541              : !> \param gxc_tau ...
     542              : ! **************************************************************************************************
     543          264 :    SUBROUTINE qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
     544              : 
     545              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     546              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
     547              : 
     548              :       INTEGER                                            :: ispin
     549              :       TYPE(pw_env_type), POINTER                         :: pw_env
     550              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     551              : 
     552          264 :       CALL get_ks_env(ks_env, pw_env=pw_env)
     553          264 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     554              : 
     555          264 :       IF (ASSOCIATED(fxc_rho)) THEN
     556          564 :          DO ispin = 1, SIZE(fxc_rho)
     557          564 :             CALL auxbas_pw_pool%give_back_pw(fxc_rho(ispin))
     558              :          END DO
     559          264 :          DEALLOCATE (fxc_rho)
     560              :       END IF
     561          264 :       IF (ASSOCIATED(fxc_tau)) THEN
     562            0 :          DO ispin = 1, SIZE(fxc_tau)
     563            0 :             CALL auxbas_pw_pool%give_back_pw(fxc_tau(ispin))
     564              :          END DO
     565            0 :          DEALLOCATE (fxc_tau)
     566              :       END IF
     567          264 :       IF (ASSOCIATED(gxc_rho)) THEN
     568          564 :          DO ispin = 1, SIZE(gxc_rho)
     569          564 :             CALL auxbas_pw_pool%give_back_pw(gxc_rho(ispin))
     570              :          END DO
     571          264 :          DEALLOCATE (gxc_rho)
     572              :       END IF
     573          264 :       IF (ASSOCIATED(gxc_tau)) THEN
     574            0 :          DO ispin = 1, SIZE(gxc_tau)
     575            0 :             CALL auxbas_pw_pool%give_back_pw(gxc_tau(ispin))
     576              :          END DO
     577            0 :          DEALLOCATE (gxc_tau)
     578              :       END IF
     579              : 
     580          264 :    END SUBROUTINE qs_fgxc_release
     581              : 
     582              : END MODULE qs_fxc
        

Generated by: LCOV version 2.0-1