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 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      649675 :    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      649675 :                  (xc_deriv_id == xc_deriv_pw)
      64      649675 :    END FUNCTION
      65             : 
      66             : ! **************************************************************************************************
      67             : !> \brief ...
      68             : !> \param pw_in ...
      69             : !> \param pw_out ...
      70             : !> \param xc_smooth_id ...
      71             : ! **************************************************************************************************
      72      133141 :    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      133141 :       CALL timeset(routineN, handle)
      82             : 
      83      265312 :       SELECT CASE (xc_smooth_id)
      84             :       CASE (xc_rho_no_smooth)
      85      132171 :          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         218 :          CALL pw_zero(pw_out)
     103             :          CALL pw_nn_smear_r(pw_in=pw_in, &
     104             :                             pw_out=pw_out, &
     105         218 :                             coeffs=nn50_coeffs)
     106             :       CASE default
     107      133141 :          CPABORT("Unsupported smoothing")
     108             :       END SELECT
     109             : 
     110      133141 :       CALL timestop(handle)
     111             : 
     112      133141 :    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      103101 :    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      412404 :       DO idir = 1, 3
     131      309303 :          CALL pw_zero(gradient(idir))
     132      412404 :          CALL xc_pw_derive(pw_r, tmp_g, gradient(idir), idir, xc_deriv_method_id, pw_g=pw_g)
     133             :       END DO
     134             : 
     135      103101 :    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        2666 :       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        2666 :          CALL timeset(routineN, handle)
     160             : 
     161        5332 :          SELECT CASE (deriv_method_id)
     162             :          CASE (xc_deriv_pw)
     163             : 
     164        2666 :             IF (PRESENT(tmp_g)) my_tmp_g = tmp_g
     165             : 
     166        2666 :             owns_tmp_g = .FALSE.
     167        2666 :             IF (.NOT. ASSOCIATED(my_tmp_g%pw_grid)) THEN
     168        1198 :                CALL pw_pool%create_pw(my_tmp_g)
     169        1198 :                owns_tmp_g = .TRUE.
     170             :             END IF
     171        2666 :             CALL pw_zero(my_tmp_g)
     172        2666 :             CALL pw_transfer(pw, my_tmp_g)
     173             : 
     174        2666 :             CALL pw_laplace(my_tmp_g)
     175             : 
     176        2666 :             IF (PRESENT(pw_out)) THEN
     177        1468 :                CALL pw_transfer(my_tmp_g, pw_out)
     178             :             ELSE
     179        1198 :                CALL pw_transfer(my_tmp_g, pw)
     180             :             END IF
     181        2666 :             IF (owns_tmp_g) THEN
     182        1198 :                CALL pw_pool%give_back_pw(my_tmp_g)
     183             :             END IF
     184             :          CASE default
     185        2666 :             CPABORT("Unsupported derivative method")
     186             :          END SELECT
     187             : 
     188        2666 :          CALL timestop(handle)
     189             : 
     190        2666 :       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       87067 :    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       87067 :       CALL timeset(routineN, handle)
     212             : 
     213             :       ! partial integration
     214       87067 :       IF (xc_deriv_method_id /= xc_deriv_pw) THEN
     215        8212 :          CALL pw_spline_scale_deriv(pw_to_deriv, transpose=.TRUE.)
     216             :       END IF
     217             : 
     218       87067 :       IF (ASSOCIATED(vxc_g%pw_grid)) CALL pw_zero(vxc_g)
     219             : 
     220      348268 :       DO idir = 1, 3
     221      261201 :          CALL xc_pw_derive(pw_to_deriv(idir), tmp_g, vxc_r, idir, xc_deriv_method_id, copy_to_vxcr=.FALSE.)
     222      348268 :          IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(vxc_g%pw_grid)) CALL pw_axpy(tmp_g, vxc_g)
     223             :       END DO
     224             : 
     225       87067 :       IF (ASSOCIATED(vxc_g%pw_grid)) THEN
     226       84637 :          CALL pw_transfer(vxc_g, pw_to_deriv(1))
     227       84637 :          CALL pw_axpy(pw_to_deriv(1), vxc_r)
     228             :       END IF
     229             : 
     230       87067 :       CALL timestop(handle)
     231             : 
     232       87067 :    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      570504 :       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      570504 :          CALL timeset(routineN, handle)
     263             : 
     264      570504 :          my_copy_to_vxcr = .TRUE.
     265      570504 :          IF (PRESENT(copy_to_vxcr)) my_copy_to_vxcr = copy_to_vxcr
     266             : 
     267      570504 :          IF (xc_requires_tmp_g(xc_deriv_method_id)) THEN
     268             : 
     269      555102 :             IF (PRESENT(pw_g)) THEN
     270      301191 :                IF (ASSOCIATED(pw_g%pw_grid)) THEN
     271      301191 :                   CALL pw_copy(pw_g, tmp_g)
     272             :                ELSE
     273           0 :                   CALL pw_transfer(pw, tmp_g)
     274             :                END IF
     275             :             ELSE
     276      253911 :                CALL pw_transfer(pw, tmp_g)
     277             :             END IF
     278             : 
     279     1071276 :             SELECT CASE (xc_deriv_method_id)
     280             :             CASE (xc_deriv_pw)
     281      516174 :                CALL pw_derive(tmp_g, nd(:, idir))
     282             :             CASE (xc_deriv_spline2)
     283       38532 :                CALL pw_spline2_interpolate_values_g(tmp_g)
     284       38532 :                CALL pw_spline2_deriv_g(tmp_g, idir=idir)
     285             :             CASE (xc_deriv_spline3)
     286         396 :                CALL pw_spline3_interpolate_values_g(tmp_g)
     287         396 :                CALL pw_spline3_deriv_g(tmp_g, idir=idir)
     288             :             CASE default
     289      555102 :                CPABORT("Unsupported deriv method")
     290             :             END SELECT
     291             : 
     292      555102 :             IF (my_copy_to_vxcr) CALL pw_transfer(tmp_g, vxc_r)
     293             :          ELSE
     294             :             #:if kind=="r3d_rs"
     295       26802 :                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       11400 :                                      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        1014 :                                      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        1230 :                                      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        1758 :                                      idir=idir)
     312             :                CASE default
     313       15402 :                   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      570504 :          CALL timestop(handle)
     342             : 
     343      570504 :       END SUBROUTINE xc_pw_derive_${kind}$
     344             :    #:endfor
     345             : 
     346             : END MODULE xc_util