LCOV - code coverage report
Current view: top level - src/xc - xc_util.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 85.4 % 89 76
Test Date: 2025-07-25 12:55:17 Functions: 87.5 % 8 7

            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 contains utility functions for the xc package
      10              : !> \par History
      11              : !>      03.2022 created [F. Stein]
      12              : !> \author Frederick Stein
      13              : ! **************************************************************************************************
      14              : MODULE xc_util
      15              :    USE pw_methods, ONLY: pw_axpy, &
      16              :                          pw_copy, &
      17              :                          pw_derive, &
      18              :                          pw_laplace, &
      19              :                          pw_transfer, &
      20              :                          pw_zero
      21              :    USE pw_pool_types, ONLY: pw_pool_type
      22              :    USE pw_spline_utils, ONLY: &
      23              :       nn10_coeffs, nn10_deriv_coeffs, nn50_coeffs, nn50_deriv_coeffs, pw_nn_deriv_r, &
      24              :       pw_nn_smear_r, pw_spline2_deriv_g, pw_spline2_interpolate_values_g, pw_spline3_deriv_g, &
      25              :       pw_spline3_interpolate_values_g, pw_spline_scale_deriv, spline2_coeffs, &
      26              :       spline2_deriv_coeffs, spline3_coeffs, spline3_deriv_coeffs
      27              :    USE pw_types, ONLY: &
      28              :       pw_c1d_gs_type, pw_r3d_rs_type
      29              :    USE xc_input_constants, ONLY: &
      30              :       xc_deriv_nn10_smooth, xc_deriv_nn50_smooth, xc_deriv_pw, xc_deriv_spline2, &
      31              :       xc_deriv_spline2_smooth, xc_deriv_spline3, xc_deriv_spline3_smooth, xc_rho_nn10, &
      32              :       xc_rho_nn50, xc_rho_no_smooth, xc_rho_spline2_smooth, xc_rho_spline3_smooth
      33              : #include "../base/base_uses.f90"
      34              : 
      35              :    IMPLICIT NONE
      36              : 
      37              :    PRIVATE
      38              : 
      39              :    PUBLIC :: xc_pw_smooth, xc_pw_laplace, xc_pw_divergence, xc_pw_derive, xc_requires_tmp_g, xc_pw_gradient
      40              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_util'
      41              : 
      42              :    INTERFACE xc_pw_derive
      43              :       MODULE PROCEDURE xc_pw_derive_r3d_rs, xc_pw_derive_c1d_gs
      44              :    END INTERFACE
      45              : 
      46              :    INTERFACE xc_pw_laplace
      47              :       MODULE PROCEDURE xc_pw_laplace_r3d_rs, xc_pw_laplace_c1d_gs
      48              :    END INTERFACE
      49              : 
      50              : CONTAINS
      51              : 
      52              : ! **************************************************************************************************
      53              : !> \brief ...
      54              : !> \param xc_deriv_id ...
      55              : !> \return ...
      56              : ! **************************************************************************************************
      57       661189 :    ELEMENTAL FUNCTION xc_requires_tmp_g(xc_deriv_id) RESULT(requires)
      58              :       INTEGER, INTENT(IN)                                :: xc_deriv_id
      59              :       LOGICAL                                            :: requires
      60              : 
      61              :       requires = (xc_deriv_id == xc_deriv_spline2) .OR. &
      62              :                  (xc_deriv_id == xc_deriv_spline3) .OR. &
      63       661189 :                  (xc_deriv_id == xc_deriv_pw)
      64       661189 :    END FUNCTION
      65              : 
      66              : ! **************************************************************************************************
      67              : !> \brief ...
      68              : !> \param pw_in ...
      69              : !> \param pw_out ...
      70              : !> \param xc_smooth_id ...
      71              : ! **************************************************************************************************
      72       133908 :    SUBROUTINE xc_pw_smooth(pw_in, pw_out, xc_smooth_id)
      73              :       TYPE(pw_r3d_rs_type), INTENT(IN)                          :: pw_in
      74              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                       :: pw_out
      75              :       INTEGER, INTENT(IN)                                :: xc_smooth_id
      76              : 
      77              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xc_pw_smooth'
      78              : 
      79              :       INTEGER                                            :: handle
      80              : 
      81       133908 :       CALL timeset(routineN, handle)
      82              : 
      83       266778 :       SELECT CASE (xc_smooth_id)
      84              :       CASE (xc_rho_no_smooth)
      85       132870 :          CALL pw_copy(pw_in, pw_out)
      86              :       CASE (xc_rho_spline2_smooth)
      87            0 :          CALL pw_zero(pw_out)
      88              :          CALL pw_nn_smear_r(pw_in=pw_in, &
      89              :                             pw_out=pw_out, &
      90            0 :                             coeffs=spline2_coeffs)
      91              :       CASE (xc_rho_spline3_smooth)
      92            0 :          CALL pw_zero(pw_out)
      93              :          CALL pw_nn_smear_r(pw_in=pw_in, &
      94              :                             pw_out=pw_out, &
      95            0 :                             coeffs=spline3_coeffs)
      96              :       CASE (xc_rho_nn10)
      97          752 :          CALL pw_zero(pw_out)
      98              :          CALL pw_nn_smear_r(pw_in=pw_in, &
      99              :                             pw_out=pw_out, &
     100          752 :                             coeffs=nn10_coeffs)
     101              :       CASE (xc_rho_nn50)
     102          286 :          CALL pw_zero(pw_out)
     103              :          CALL pw_nn_smear_r(pw_in=pw_in, &
     104              :                             pw_out=pw_out, &
     105          286 :                             coeffs=nn50_coeffs)
     106              :       CASE default
     107       133908 :          CPABORT("Unsupported smoothing")
     108              :       END SELECT
     109              : 
     110       133908 :       CALL timestop(handle)
     111              : 
     112       133908 :    END SUBROUTINE xc_pw_smooth
     113              : 
     114              : ! **************************************************************************************************
     115              : !> \brief ...
     116              : !> \param pw_r ...
     117              : !> \param pw_g ...
     118              : !> \param tmp_g ...
     119              : !> \param gradient ...
     120              : !> \param xc_deriv_method_id ...
     121              : ! **************************************************************************************************
     122       105721 :    SUBROUTINE xc_pw_gradient(pw_r, pw_g, tmp_g, gradient, xc_deriv_method_id)
     123              :       TYPE(pw_r3d_rs_type), INTENT(IN)                          :: pw_r
     124              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                   :: pw_g, tmp_g
     125              :       TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT)         :: gradient
     126              :       INTEGER, INTENT(IN)                                :: xc_deriv_method_id
     127              : 
     128              :       INTEGER                                            :: idir
     129              : 
     130       422884 :       DO idir = 1, 3
     131       317163 :          CALL pw_zero(gradient(idir))
     132       422884 :          CALL xc_pw_derive(pw_r, tmp_g, gradient(idir), idir, xc_deriv_method_id, pw_g=pw_g)
     133              :       END DO
     134              : 
     135       105721 :    END SUBROUTINE xc_pw_gradient
     136              : 
     137              :    #:for kind in ["r3d_rs", "c1d_gs"]
     138              : ! **************************************************************************************************
     139              : !> \brief Calculates the Laplacian of pw
     140              : !> \param pw on input: pw of which the Laplacian shall be calculated, on output if pw_out is absent: Laplacian of input
     141              : !> \param pw_pool ...
     142              : !> \param deriv_method_id ...
     143              : !> \param pw_out if present, save the Laplacian of pw here
     144              : !> \param tmp_g scratch grid in reciprocal space, used instead of the internal grid if given explicitly to save memory
     145              : ! **************************************************************************************************
     146         2982 :       SUBROUTINE xc_pw_laplace_${kind}$ (pw, pw_pool, deriv_method_id, pw_out, tmp_g)
     147              :          TYPE(pw_${kind}$_type), INTENT(INOUT)                       :: pw
     148              :          TYPE(pw_pool_type), INTENT(IN), POINTER            :: pw_pool
     149              :          INTEGER, INTENT(IN)                                :: deriv_method_id
     150              :          TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL             :: pw_out
     151              :          TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL            :: tmp_g
     152              : 
     153              :          CHARACTER(len=*), PARAMETER                        :: routineN = 'xc_pw_laplace'
     154              : 
     155              :          INTEGER                                            :: handle
     156              :          LOGICAL                                            :: owns_tmp_g
     157              :          TYPE(pw_c1d_gs_type)                                  :: my_tmp_g
     158              : 
     159         2982 :          CALL timeset(routineN, handle)
     160              : 
     161         5964 :          SELECT CASE (deriv_method_id)
     162              :          CASE (xc_deriv_pw)
     163              : 
     164         2982 :             IF (PRESENT(tmp_g)) my_tmp_g = tmp_g
     165              : 
     166         2982 :             owns_tmp_g = .FALSE.
     167         2982 :             IF (.NOT. ASSOCIATED(my_tmp_g%pw_grid)) THEN
     168         1304 :                CALL pw_pool%create_pw(my_tmp_g)
     169         1304 :                owns_tmp_g = .TRUE.
     170              :             END IF
     171         2982 :             CALL pw_zero(my_tmp_g)
     172         2982 :             CALL pw_transfer(pw, my_tmp_g)
     173              : 
     174         2982 :             CALL pw_laplace(my_tmp_g)
     175              : 
     176         2982 :             IF (PRESENT(pw_out)) THEN
     177         1678 :                CALL pw_transfer(my_tmp_g, pw_out)
     178              :             ELSE
     179         1304 :                CALL pw_transfer(my_tmp_g, pw)
     180              :             END IF
     181         2982 :             IF (owns_tmp_g) THEN
     182         1304 :                CALL pw_pool%give_back_pw(my_tmp_g)
     183              :             END IF
     184              :          CASE default
     185         2982 :             CPABORT("Unsupported derivative method")
     186              :          END SELECT
     187              : 
     188         2982 :          CALL timestop(handle)
     189              : 
     190         2982 :       END SUBROUTINE xc_pw_laplace_${kind}$
     191              :    #:endfor
     192              : 
     193              : ! **************************************************************************************************
     194              : !> \brief Calculates the divergence of pw_to_deriv
     195              : !> \param xc_deriv_method_id ...
     196              : !> \param pw_to_deriv ...
     197              : !> \param tmp_g ...
     198              : !> \param vxc_g ...
     199              : !> \param vxc_r ...
     200              : ! **************************************************************************************************
     201        88005 :    SUBROUTINE xc_pw_divergence(xc_deriv_method_id, pw_to_deriv, tmp_g, vxc_g, vxc_r)
     202              :       INTEGER, INTENT(IN)                                :: xc_deriv_method_id
     203              :       TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT)         :: pw_to_deriv
     204              :       TYPE(pw_c1d_gs_type), INTENT(INOUT)                   :: tmp_g, vxc_g
     205              :       TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vxc_r
     206              : 
     207              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xc_pw_divergence'
     208              : 
     209              :       INTEGER                                            :: handle, idir
     210              : 
     211        88005 :       CALL timeset(routineN, handle)
     212              : 
     213              :       ! partial integration
     214        88005 :       IF (xc_deriv_method_id /= xc_deriv_pw) THEN
     215         7556 :          CALL pw_spline_scale_deriv(pw_to_deriv, transpose=.TRUE.)
     216              :       END IF
     217              : 
     218        88005 :       IF (ASSOCIATED(vxc_g%pw_grid)) CALL pw_zero(vxc_g)
     219              : 
     220       352020 :       DO idir = 1, 3
     221       264015 :          CALL xc_pw_derive(pw_to_deriv(idir), tmp_g, vxc_r, idir, xc_deriv_method_id, copy_to_vxcr=.FALSE.)
     222       352020 :          IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(vxc_g%pw_grid)) CALL pw_axpy(tmp_g, vxc_g)
     223              :       END DO
     224              : 
     225        88005 :       IF (ASSOCIATED(vxc_g%pw_grid)) THEN
     226        86225 :          CALL pw_transfer(vxc_g, pw_to_deriv(1))
     227        86225 :          CALL pw_axpy(pw_to_deriv(1), vxc_r)
     228              :       END IF
     229              : 
     230        88005 :       CALL timestop(handle)
     231              : 
     232        88005 :    END SUBROUTINE xc_pw_divergence
     233              : 
     234              :    #:for kind in ["r3d_rs", "c1d_gs"]
     235              : ! **************************************************************************************************
     236              : !> \brief Calculates the derivative of a function on a planewave grid in a given direction
     237              : !> \param pw function to derive
     238              : !> \param tmp_g temporary grid in reciprocal space, only required if derivative method is pw or spline
     239              : !> \param vxc_r if tmp_g is not required, add derivative here
     240              : !> \param idir direction of derivative
     241              : !> \param xc_deriv_method_id ...
     242              : !> \param copy_to_vxcr ...
     243              : !> \param pw_g ...
     244              : ! **************************************************************************************************
     245       581178 :       SUBROUTINE xc_pw_derive_${kind}$ (pw, tmp_g, vxc_r, idir, xc_deriv_method_id, copy_to_vxcr, pw_g)
     246              :          TYPE(pw_${kind}$_type), INTENT(IN)                          :: pw
     247              :          TYPE(pw_c1d_gs_type), INTENT(INOUT)                   :: tmp_g
     248              :          TYPE(pw_r3d_rs_type), INTENT(INOUT)                       :: vxc_r
     249              :          INTEGER, INTENT(IN)                                :: idir, xc_deriv_method_id
     250              :          LOGICAL, INTENT(IN), OPTIONAL                      :: copy_to_vxcr
     251              :          TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL            :: pw_g
     252              : 
     253              :          CHARACTER(len=*), PARAMETER                        :: routineN = 'xc_pw_derive'
     254              :          INTEGER, DIMENSION(3, 3), PARAMETER :: nd = RESHAPE((/1, 0, 0, 0, 1, 0, 0, 0, 1/), (/3, 3/))
     255              : 
     256              :          INTEGER                                            :: handle
     257              :          LOGICAL                                            :: my_copy_to_vxcr
     258              :          #:if kind=="c1d_gs"
     259              :             TYPE(pw_r3d_rs_type) :: tmp_r
     260              :          #:endif
     261              : 
     262       581178 :          CALL timeset(routineN, handle)
     263              : 
     264       581178 :          my_copy_to_vxcr = .TRUE.
     265       581178 :          IF (PRESENT(copy_to_vxcr)) my_copy_to_vxcr = copy_to_vxcr
     266              : 
     267       581178 :          IF (xc_requires_tmp_g(xc_deriv_method_id)) THEN
     268              : 
     269       569700 :             IF (PRESENT(pw_g)) THEN
     270       311025 :                IF (ASSOCIATED(pw_g%pw_grid)) THEN
     271       311025 :                   CALL pw_copy(pw_g, tmp_g)
     272              :                ELSE
     273            0 :                   CALL pw_transfer(pw, tmp_g)
     274              :                END IF
     275              :             ELSE
     276       258675 :                CALL pw_transfer(pw, tmp_g)
     277              :             END IF
     278              : 
     279      1100508 :             SELECT CASE (xc_deriv_method_id)
     280              :             CASE (xc_deriv_pw)
     281       530808 :                CALL pw_derive(tmp_g, nd(:, idir))
     282              :             CASE (xc_deriv_spline2)
     283        38550 :                CALL pw_spline2_interpolate_values_g(tmp_g)
     284        38550 :                CALL pw_spline2_deriv_g(tmp_g, idir=idir)
     285              :             CASE (xc_deriv_spline3)
     286          342 :                CALL pw_spline3_interpolate_values_g(tmp_g)
     287          342 :                CALL pw_spline3_deriv_g(tmp_g, idir=idir)
     288              :             CASE default
     289       569700 :                CPABORT("Unsupported deriv method")
     290              :             END SELECT
     291              : 
     292       569700 :             IF (my_copy_to_vxcr) CALL pw_transfer(tmp_g, vxc_r)
     293              :          ELSE
     294              :             #:if kind=="r3d_rs"
     295        18666 :                SELECT CASE (xc_deriv_method_id)
     296              :                CASE (xc_deriv_spline2_smooth)
     297              :                   CALL pw_nn_deriv_r(pw_in=pw, &
     298              :                                      pw_out=vxc_r, coeffs=spline2_deriv_coeffs, &
     299         7188 :                                      idir=idir)
     300              :                CASE (xc_deriv_spline3_smooth)
     301              :                   CALL pw_nn_deriv_r(pw_in=pw, &
     302              :                                      pw_out=vxc_r, coeffs=spline3_deriv_coeffs, &
     303         1098 :                                      idir=idir)
     304              :                CASE (xc_deriv_nn10_smooth)
     305              :                   CALL pw_nn_deriv_r(pw_in=pw, &
     306              :                                      pw_out=vxc_r, coeffs=nn10_deriv_coeffs, &
     307         1284 :                                      idir=idir)
     308              :                CASE (xc_deriv_nn50_smooth)
     309              :                   CALL pw_nn_deriv_r(pw_in=pw, &
     310              :                                      pw_out=vxc_r, coeffs=nn50_deriv_coeffs, &
     311         1908 :                                      idir=idir)
     312              :                CASE default
     313        11478 :                   CPABORT("Unsupported derivative method")
     314              :                END SELECT
     315              :             #:else
     316            0 :                CALL tmp_r%create(pw%pw_grid)
     317            0 :                SELECT CASE (xc_deriv_method_id)
     318              :                CASE (xc_deriv_spline2_smooth)
     319              :                   CALL pw_nn_deriv_r(pw_in=tmp_r, &
     320              :                                      pw_out=vxc_r, coeffs=spline2_deriv_coeffs, &
     321            0 :                                      idir=idir)
     322              :                CASE (xc_deriv_spline3_smooth)
     323              :                   CALL pw_nn_deriv_r(pw_in=tmp_r, &
     324              :                                      pw_out=vxc_r, coeffs=spline3_deriv_coeffs, &
     325            0 :                                      idir=idir)
     326              :                CASE (xc_deriv_nn10_smooth)
     327              :                   CALL pw_nn_deriv_r(pw_in=tmp_r, &
     328              :                                      pw_out=vxc_r, coeffs=nn10_deriv_coeffs, &
     329            0 :                                      idir=idir)
     330              :                CASE (xc_deriv_nn50_smooth)
     331              :                   CALL pw_nn_deriv_r(pw_in=tmp_r, &
     332              :                                      pw_out=vxc_r, coeffs=nn50_deriv_coeffs, &
     333            0 :                                      idir=idir)
     334              :                CASE default
     335            0 :                   CPABORT("Unsupported derivative method")
     336              :                END SELECT
     337            0 :                CALL tmp_r%release()
     338              :             #:endif
     339              :          END IF
     340              : 
     341       581178 :          CALL timestop(handle)
     342              : 
     343       581178 :       END SUBROUTINE xc_pw_derive_${kind}$
     344              :    #:endfor
     345              : 
     346              : END MODULE xc_util
        

Generated by: LCOV version 2.0-1