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 Setting up the Spline coefficients used to Interpolate the G-Term
      10             : !>      in Ewald sums
      11             : !> \par History
      12             : !>      12.2005 created [tlaino]
      13             : !> \author Teodoro Laino
      14             : ! **************************************************************************************************
      15             : MODULE ewald_spline_util
      16             :    USE cell_methods,                    ONLY: cell_create
      17             :    USE cell_types,                      ONLY: cell_release,&
      18             :                                               cell_type
      19             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      20             :                                               cp_logger_type
      21             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      22             :                                               cp_print_key_unit_nr
      23             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      24             :                                               section_vals_type,&
      25             :                                               section_vals_val_get
      26             :    USE kinds,                           ONLY: dp
      27             :    USE message_passing,                 ONLY: mp_para_env_type
      28             :    USE pw_grid_types,                   ONLY: HALFSPACE,&
      29             :                                               pw_grid_type
      30             :    USE pw_grids,                        ONLY: pw_grid_create,&
      31             :                                               pw_grid_setup
      32             :    USE pw_methods,                      ONLY: pw_zero
      33             :    USE pw_pool_types,                   ONLY: pw_pool_create,&
      34             :                                               pw_pool_type
      35             :    USE pw_spline_utils,                 ONLY: &
      36             :         Eval_Interp_Spl3_pbc, Eval_d_Interp_Spl3_pbc, find_coeffs, pw_spline_do_precond, &
      37             :         pw_spline_precond_create, pw_spline_precond_release, pw_spline_precond_set_kind, &
      38             :         pw_spline_precond_type, spl3_pbc
      39             :    USE pw_types,                        ONLY: pw_r3d_rs_type
      40             : !NB parallelization
      41             : #include "./base/base_uses.f90"
      42             : 
      43             :    IMPLICIT NONE
      44             :    PRIVATE
      45             : 
      46             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      47             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_spline_util'
      48             :    PUBLIC :: Setup_Ewald_Spline
      49             : 
      50             : CONTAINS
      51             : 
      52             : ! **************************************************************************************************
      53             : !> \brief Setup of the G-space Ewald Term Spline Coefficients
      54             : !> \param pw_grid ...
      55             : !> \param pw_pool ...
      56             : !> \param coeff ...
      57             : !> \param LG ...
      58             : !> \param gx ...
      59             : !> \param gy ...
      60             : !> \param gz ...
      61             : !> \param hmat ...
      62             : !> \param npts ...
      63             : !> \param param_section ...
      64             : !> \param tag ...
      65             : !> \param print_section ...
      66             : !> \param para_env ...
      67             : !> \par History
      68             : !>      12.2005 created [tlaino]
      69             : !> \author Teodoro Laino
      70             : ! **************************************************************************************************
      71         102 :    SUBROUTINE Setup_Ewald_Spline(pw_grid, pw_pool, coeff, LG, gx, gy, gz, hmat, npts, &
      72             :                                  param_section, tag, print_section, para_env)
      73             :       TYPE(pw_grid_type), POINTER                        :: pw_grid
      74             :       TYPE(pw_pool_type), POINTER                        :: pw_pool
      75             :       TYPE(pw_r3d_rs_type), POINTER                      :: coeff
      76             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: LG, gx, gy, gz
      77             :       REAL(KIND=dp), INTENT(IN)                          :: hmat(3, 3)
      78             :       INTEGER, INTENT(IN)                                :: npts(3)
      79             :       TYPE(section_vals_type), POINTER                   :: param_section
      80             :       CHARACTER(LEN=*), INTENT(IN)                       :: tag
      81             :       TYPE(section_vals_type), POINTER                   :: print_section
      82             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      83             : 
      84             :       INTEGER                                            :: bo(2, 3), iounit
      85             :       TYPE(cell_type), POINTER                           :: cell
      86             :       TYPE(cp_logger_type), POINTER                      :: logger
      87             :       TYPE(pw_r3d_rs_type)                               :: pw
      88             : 
      89             : !
      90             : ! Setting Up Fit Procedure
      91             : !
      92             : 
      93           0 :       CPASSERT(.NOT. ASSOCIATED(pw_grid))
      94         102 :       CPASSERT(.NOT. ASSOCIATED(pw_pool))
      95         102 :       CPASSERT(.NOT. ASSOCIATED(coeff))
      96         102 :       NULLIFY (cell)
      97             : 
      98         102 :       CALL cell_create(cell, hmat=hmat, periodic=(/1, 1, 1/))
      99         102 :       CALL pw_grid_create(pw_grid, para_env, local=.TRUE.)
     100         102 :       logger => cp_get_default_logger()
     101             :       iounit = cp_print_key_unit_nr(logger, print_section, "", &
     102         102 :                                     extension=".Log")
     103         408 :       bo(1, 1:3) = 0
     104         408 :       bo(2, 1:3) = npts(1:3) - 1
     105         102 :       CALL pw_grid_setup(cell%hmat, pw_grid, grid_span=HALFSPACE, bounds=bo, iounit=iounit)
     106             : 
     107             :       CALL cp_print_key_finished_output(iounit, logger, print_section, &
     108         102 :                                         "")
     109             :       ! pw_pool initialized
     110         102 :       CALL pw_pool_create(pw_pool, pw_grid=pw_grid)
     111         102 :       ALLOCATE (coeff)
     112         102 :       CALL pw_pool%create_pw(pw)
     113         102 :       CALL pw_pool%create_pw(coeff)
     114             :       ! Evaluate function on grid
     115             :       CALL eval_pw_TabLR(pw, pw_pool, coeff, Lg, gx, gy, gz, hmat_mm=hmat, &
     116         102 :                          param_section=param_section, tag=tag)
     117         102 :       CALL pw_pool%give_back_pw(pw)
     118         102 :       CALL cell_release(cell)
     119             : 
     120         102 :    END SUBROUTINE Setup_Ewald_Spline
     121             : 
     122             : ! **************************************************************************************************
     123             : !> \brief Evaluates the function G-Term in reciprocal space on the grid
     124             : !>      and find the coefficients of the Splines
     125             : !> \param grid ...
     126             : !> \param pw_pool ...
     127             : !> \param TabLR ...
     128             : !> \param Lg ...
     129             : !> \param gx ...
     130             : !> \param gy ...
     131             : !> \param gz ...
     132             : !> \param hmat_mm ...
     133             : !> \param param_section ...
     134             : !> \param tag ...
     135             : !> \par History
     136             : !>      12.2005 created [tlaino]
     137             : !> \author Teodoro Laino
     138             : ! **************************************************************************************************
     139         102 :    SUBROUTINE eval_pw_TabLR(grid, pw_pool, TabLR, Lg, gx, gy, gz, hmat_mm, &
     140             :                             param_section, tag)
     141             :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: grid
     142             :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     143             :       TYPE(pw_r3d_rs_type), POINTER                      :: TabLR
     144             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: Lg, gx, gy, gz
     145             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat_mm
     146             :       TYPE(section_vals_type), POINTER                   :: param_section
     147             :       CHARACTER(LEN=*), INTENT(IN)                       :: tag
     148             : 
     149             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'eval_pw_TabLR'
     150             : 
     151             :       INTEGER :: act_nx, act_ny, aint_precond, handle, i, iii, is, j, js, k, ks, Lg_loc_max, &
     152             :          Lg_loc_min, max_iter, my_i, my_j, my_k, n1, n2, n3, n_extra, NLg_loc, nxlim, nylim, &
     153             :          nzlim, precond_kind
     154             :       INTEGER, DIMENSION(2, 3)                           :: gbo
     155             :       LOGICAL                                            :: success
     156             :       REAL(KIND=dp)                                      :: dr1, dr2, dr3, eps_r, eps_x, xs1, xs2, &
     157             :                                                             xs3
     158         102 :       REAL(KIND=dp), ALLOCATABLE                         :: cos_gx(:, :), cos_gy(:, :), &
     159         102 :                                                             cos_gz(:, :), lhs(:, :), rhs(:, :), &
     160         102 :                                                             sin_gx(:, :), sin_gy(:, :), &
     161         102 :                                                             sin_gz(:, :)
     162             :       TYPE(pw_spline_precond_type)                       :: precond
     163             :       TYPE(section_vals_type), POINTER                   :: interp_section
     164             : 
     165             : !NB pull expensive Cos() out of inner looop
     166             : !NB temporaries for holding stuff so that dgemm can be used
     167             : 
     168             :       EXTERNAL :: DGEMM
     169             : 
     170         102 :       CALL timeset(routineN, handle)
     171         102 :       n1 = grid%pw_grid%npts(1)
     172         102 :       n2 = grid%pw_grid%npts(2)
     173         102 :       n3 = grid%pw_grid%npts(3)
     174         102 :       dr1 = grid%pw_grid%dr(1)
     175         102 :       dr2 = grid%pw_grid%dr(2)
     176         102 :       dr3 = grid%pw_grid%dr(3)
     177        1020 :       gbo = grid%pw_grid%bounds
     178         102 :       nxlim = FLOOR(REAL(n1, KIND=dp)/2.0_dp)
     179         102 :       nylim = FLOOR(REAL(n2, KIND=dp)/2.0_dp)
     180         102 :       nzlim = FLOOR(REAL(n3, KIND=dp)/2.0_dp)
     181         102 :       is = 0
     182         102 :       js = 0
     183         102 :       ks = 0
     184         102 :       IF (2*nxlim /= n1) is = 1
     185         102 :       IF (2*nylim /= n2) js = 1
     186         102 :       IF (2*nzlim /= n3) ks = 1
     187         102 :       CALL pw_zero(grid)
     188             : 
     189             :       ! Used the full symmetry to reduce the evaluation to 1/64th
     190             :       !NB parallelization
     191         102 :       iii = 0
     192             :       !NB allocate temporaries for Cos refactoring
     193         408 :       ALLOCATE (cos_gx(SIZE(Lg), gbo(1, 1):gbo(2, 1)))
     194         408 :       ALLOCATE (sin_gx(SIZE(Lg), gbo(1, 1):gbo(2, 1)))
     195         408 :       ALLOCATE (cos_gy(SIZE(Lg), gbo(1, 2):gbo(2, 2)))
     196         408 :       ALLOCATE (sin_gy(SIZE(Lg), gbo(1, 2):gbo(2, 2)))
     197         408 :       ALLOCATE (cos_gz(SIZE(Lg), gbo(1, 3):gbo(2, 3)))
     198         408 :       ALLOCATE (sin_gz(SIZE(Lg), gbo(1, 3):gbo(2, 3)))
     199             :       !NB precalculate Cos(gx*xs1) etc for Cos refactoring
     200        5094 :       DO k = gbo(1, 3), gbo(2, 3)
     201        4992 :          my_k = k - gbo(1, 3)
     202        4992 :          xs3 = REAL(my_k, dp)*dr3
     203        4992 :          IF (k > nzlim) CYCLE
     204     1824444 :          cos_gz(1:SIZE(Lg), k) = COS(gz(1:SIZE(Lg))*xs3)
     205     1824546 :          sin_gz(1:SIZE(Lg), k) = SIN(gz(1:SIZE(Lg))*xs3)
     206             :       END DO ! k
     207             :       xs2 = 0.0_dp
     208        5118 :       DO j = gbo(1, 2), gbo(2, 2)
     209        5016 :          IF (j > nylim) CYCLE
     210     1920156 :          cos_gy(1:SIZE(Lg), j) = COS(gy(1:SIZE(Lg))*xs2)
     211     1920156 :          sin_gy(1:SIZE(Lg), j) = SIN(gy(1:SIZE(Lg))*xs2)
     212        5118 :          xs2 = xs2 + dr2
     213             :       END DO ! j
     214             :       xs1 = 0.0_dp
     215        5070 :       DO i = gbo(1, 1), gbo(2, 1)
     216        4968 :          IF (i > nxlim) CYCLE
     217     1728732 :          cos_gx(1:SIZE(Lg), i) = COS(gx(1:SIZE(Lg))*xs1)
     218     1728732 :          sin_gx(1:SIZE(Lg), i) = SIN(gx(1:SIZE(Lg))*xs1)
     219        5070 :          xs1 = xs1 + dr1
     220             :       END DO ! i
     221             : 
     222             :       !NB use DGEMM to compute sum over kg for each i, j, k
     223             :       ! number of elements per node, round down
     224         102 :       NLg_loc = SIZE(Lg)/grid%pw_grid%para%group_size
     225             :       ! number of extra elements not yet accounted for
     226         102 :       n_extra = MOD(SIZE(Lg), grid%pw_grid%para%group_size)
     227             :       ! first n_extra nodes get NLg_loc+1, remaining get NLg_loc
     228         102 :       IF (grid%pw_grid%para%my_pos < n_extra) THEN
     229           0 :          Lg_loc_min = (NLg_loc + 1)*grid%pw_grid%para%my_pos + 1
     230           0 :          Lg_loc_max = Lg_loc_min + (NLg_loc + 1) - 1
     231             :       ELSE
     232         102 :          Lg_loc_min = (NLg_loc + 1)*n_extra + NLg_loc*(grid%pw_grid%para%my_pos - n_extra) + 1
     233         102 :          Lg_loc_max = Lg_loc_min + NLg_loc - 1
     234             :       END IF
     235             :       ! shouldn't be necessary
     236         102 :       Lg_loc_max = MIN(SIZE(Lg), Lg_loc_max)
     237         102 :       NLg_loc = Lg_loc_max - Lg_loc_min + 1
     238             : 
     239         102 :       IF (NLg_loc > 0) THEN ! some work for this node
     240             : 
     241          55 :          act_nx = MIN(gbo(2, 1), nxlim) - gbo(1, 1) + 1
     242          55 :          act_ny = MIN(gbo(2, 2), nylim) - gbo(1, 2) + 1
     243             :          !NB temporaries for DGEMM use
     244         220 :          ALLOCATE (lhs(act_nx, NLg_loc))
     245         220 :          ALLOCATE (rhs(act_ny, NLg_loc))
     246             : 
     247             :          ! do cos(gx) cos(gy+gz) term
     248        2739 :          DO i = gbo(1, 1), gbo(2, 1)
     249        2684 :             IF (i > nxlim) CYCLE
     250      940029 :             lhs(i - gbo(1, 1) + 1, 1:NLg_loc) = lg(Lg_loc_min:Lg_loc_max)*cos_gx(Lg_loc_min:Lg_loc_max, i)
     251             :          END DO
     252        2751 :          DO k = gbo(1, 3), gbo(2, 3)
     253        2696 :             IF (k > nzlim) CYCLE
     254       70998 :             DO j = gbo(1, 2), gbo(2, 2)
     255       69596 :                IF (j > nylim) CYCLE
     256             :                rhs(j - gbo(1, 2) + 1, 1:NLg_loc) = cos_gy(Lg_loc_min:Lg_loc_max, j)*cos_gz(Lg_loc_min:Lg_loc_max, k) - &
     257    23579318 :                                                    sin_gy(Lg_loc_min:Lg_loc_max, j)*sin_gz(Lg_loc_min:Lg_loc_max, k)
     258             :             END DO
     259             :             CALL DGEMM('N', 'T', act_nx, act_ny, NLg_loc, 1.0D0, lhs(1, 1), act_nx, rhs(1, 1), act_ny, 0.0D0, &
     260        2751 :                        grid%array(gbo(1, 1), gbo(1, 2), k), SIZE(grid%array, 1))
     261             :          END DO
     262             : 
     263             :          ! do sin(gx) sin(gy+gz) term
     264        2739 :          DO i = gbo(1, 1), gbo(2, 1)
     265        2684 :             IF (i > nxlim) CYCLE
     266      940029 :             lhs(i - gbo(1, 1) + 1, 1:NLg_loc) = -lg(Lg_loc_min:Lg_loc_max)*sin_gx(Lg_loc_min:Lg_loc_max, i)
     267             :          END DO
     268        2751 :          DO k = gbo(1, 3), gbo(2, 3)
     269        2696 :             IF (k > nzlim) CYCLE
     270       70998 :             DO j = gbo(1, 2), gbo(2, 2)
     271       69596 :                IF (j > nylim) CYCLE
     272             :                rhs(j - gbo(1, 2) + 1, 1:NLg_loc) = cos_gy(Lg_loc_min:Lg_loc_max, j)*sin_gz(Lg_loc_min:Lg_loc_max, k) + &
     273    23579318 :                                                    sin_gy(Lg_loc_min:Lg_loc_max, j)*cos_gz(Lg_loc_min:Lg_loc_max, k)
     274             :             END DO
     275             :             CALL DGEMM('N', 'T', act_nx, act_ny, NLg_loc, 1.0D0, lhs(1, 1), act_nx, rhs(1, 1), act_ny, 1.0D0, &
     276        2751 :                        grid%array(gbo(1, 1), gbo(1, 2), k), SIZE(grid%array, 1))
     277             :          END DO
     278             : 
     279             :          !NB deallocate temporaries for DGEMM use
     280          55 :          DEALLOCATE (lhs)
     281          55 :          DEALLOCATE (rhs)
     282             :          !NB deallocate temporaries for Cos refactoring
     283          55 :          DEALLOCATE (cos_gx)
     284          55 :          DEALLOCATE (sin_gx)
     285          55 :          DEALLOCATE (cos_gy)
     286          55 :          DEALLOCATE (sin_gy)
     287          55 :          DEALLOCATE (cos_gz)
     288          55 :          DEALLOCATE (sin_gz)
     289             :          !NB parallelization
     290             :       ELSE ! no work for this node, just zero contribution
     291      826181 :          grid%array(gbo(1, 1):nxlim, gbo(1, 2):nylim, gbo(1, 3):nzlim) = 0.0_dp
     292             :       END IF ! NLg_loc > 0
     293             : 
     294     3597086 :       CALL grid%pw_grid%para%group%sum(grid%array(gbo(1, 1):nxlim, gbo(1, 2):nylim, gbo(1, 3):nzlim))
     295             : 
     296        5094 :       Fake_LoopOnGrid: DO k = gbo(1, 3), gbo(2, 3)
     297        4992 :          my_k = k
     298        4992 :          IF (k > nzlim) my_k = nzlim - ABS(nzlim - k) + ks
     299      252762 :          DO j = gbo(1, 2), gbo(2, 2)
     300      247668 :             my_j = j
     301      247668 :             IF (j > nylim) my_j = nylim - ABS(nylim - j) + js
     302    12548016 :             DO i = gbo(1, 1), gbo(2, 1)
     303    12295356 :                my_i = i
     304    12295356 :                IF (i > nxlim) my_i = nxlim - ABS(nxlim - i) + is
     305    12543024 :                grid%array(i, j, k) = grid%array(my_i, my_j, my_k)
     306             :             END DO
     307             :          END DO
     308             :       END DO Fake_LoopOnGrid
     309             :       !
     310             :       ! Solve for spline coefficients
     311             :       !
     312         102 :       interp_section => section_vals_get_subs_vals(param_section, "INTERPOLATOR")
     313         102 :       CALL section_vals_val_get(interp_section, "aint_precond", i_val=aint_precond)
     314         102 :       CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind)
     315         102 :       CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter)
     316         102 :       CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r)
     317         102 :       CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x)
     318             :       !
     319             :       ! Solve for spline coefficients
     320             :       !
     321             :       CALL pw_spline_precond_create(precond, precond_kind=aint_precond, &
     322         102 :                                     pool=pw_pool, pbc=.TRUE., transpose=.FALSE.)
     323         102 :       CALL pw_spline_do_precond(precond, grid, TabLR)
     324         102 :       CALL pw_spline_precond_set_kind(precond, precond_kind)
     325             :       success = find_coeffs(values=grid, coeffs=TabLR, &
     326             :                             linOp=spl3_pbc, preconditioner=precond, pool=pw_pool, &
     327             :                             eps_r=eps_r, eps_x=eps_x, &
     328         102 :                             max_iter=max_iter)
     329         102 :       CPASSERT(success)
     330         102 :       CALL pw_spline_precond_release(precond)
     331             :       !
     332             :       ! Check for the interpolation Spline
     333             :       !
     334             :       CALL check_spline_interp_TabLR(hmat_mm, Lg, gx, gy, gz, TabLR, param_section, &
     335         102 :                                      tag)
     336         102 :       CALL timestop(handle)
     337        1020 :    END SUBROUTINE eval_pw_TabLR
     338             : 
     339             : ! **************************************************************************************************
     340             : !> \brief Routine to check the accuracy for the Spline Interpolation
     341             : !> \param hmat_mm ...
     342             : !> \param Lg ...
     343             : !> \param gx ...
     344             : !> \param gy ...
     345             : !> \param gz ...
     346             : !> \param TabLR ...
     347             : !> \param param_section ...
     348             : !> \param tag ...
     349             : !> \par History
     350             : !>      12.2005 created [tlaino]
     351             : !> \author Teodoro Laino
     352             : ! **************************************************************************************************
     353         204 :    SUBROUTINE check_spline_interp_TabLR(hmat_mm, Lg, gx, gy, gz, TabLR, &
     354             :                                         param_section, tag)
     355             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat_mm
     356             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: Lg, gx, gy, gz
     357             :       TYPE(pw_r3d_rs_type), POINTER                      :: TabLR
     358             :       TYPE(section_vals_type), POINTER                   :: param_section
     359             :       CHARACTER(LEN=*), INTENT(IN)                       :: tag
     360             : 
     361             :       CHARACTER(len=*), PARAMETER :: routineN = 'check_spline_interp_TabLR'
     362             : 
     363             :       INTEGER                                            :: handle, i, iw, kg, npoints
     364             :       REAL(KIND=dp)                                      :: dn(3), dr1, dr2, dr3, dxTerm, dyTerm, &
     365             :                                                             dzTerm, errd, errf, Fterm, maxerrord, &
     366             :                                                             maxerrorf, Na, Nn, Term, tmp1, tmp2, &
     367             :                                                             vec(3), xs1, xs2, xs3
     368             :       TYPE(cp_logger_type), POINTER                      :: logger
     369             : 
     370         102 :       NULLIFY (logger)
     371         102 :       logger => cp_get_default_logger()
     372             :       iw = cp_print_key_unit_nr(logger, param_section, "check_spline", &
     373         102 :                                 extension="."//TRIM(tag)//"Log")
     374         102 :       CALL timeset(routineN, handle)
     375         102 :       IF (iw > 0) THEN
     376          46 :          npoints = 100
     377          46 :          errf = 0.0_dp
     378          46 :          maxerrorf = 0.0_dp
     379          46 :          errd = 0.0_dp
     380          46 :          maxerrord = 0.0_dp
     381          46 :          dr1 = hmat_mm(1, 1)/REAL(npoints, KIND=dp)
     382          46 :          dr2 = hmat_mm(2, 2)/REAL(npoints, KIND=dp)
     383          46 :          dr3 = hmat_mm(3, 3)/REAL(npoints, KIND=dp)
     384          46 :          xs1 = 0.0_dp
     385          46 :          xs2 = 0.0_dp
     386          46 :          xs3 = 0.0_dp
     387             :          WRITE (iw, '(A,T5,A15,4X,A17,T50,4X,A,5X,A,T80,A,T85,A15,4X,A17,T130,4X,A,5X,A)') &
     388          46 :             "#", "Analytical Term", "Interpolated Term", "Error", "MaxError", &
     389          92 :             "*", " Analyt Deriv  ", "Interp Deriv Mod ", "Error", "MaxError"
     390        4692 :          DO i = 1, npoints + 1
     391        4646 :             Term = 0.0_dp
     392        4646 :             dxTerm = 0.0_dp
     393        4646 :             dyTerm = 0.0_dp
     394        4646 :             dzTerm = 0.0_dp
     395             :             ! Sum over k vectors
     396     4670947 :             DO kg = 1, SIZE(Lg)
     397    18665204 :                vec = (/REAL(gx(kg), KIND=dp), REAL(gy(kg), KIND=dp), REAL(gz(kg), KIND=dp)/)
     398     4666301 :                Term = Term + lg(kg)*COS(vec(1)*xs1 + vec(2)*xs2 + vec(3)*xs3)
     399     4666301 :                dxTerm = dxTerm - lg(kg)*SIN(vec(1)*xs1 + vec(2)*xs2 + vec(3)*xs3)*vec(1)
     400     4666301 :                dyTerm = dyTerm - lg(kg)*SIN(vec(1)*xs1 + vec(2)*xs2 + vec(3)*xs3)*vec(2)
     401     4670947 :                dzTerm = dzTerm - lg(kg)*SIN(vec(1)*xs1 + vec(2)*xs2 + vec(3)*xs3)*vec(3)
     402             :             END DO
     403        4646 :             Na = SQRT(dxTerm*dxTerm + dyTerm*dyTerm + dzTerm*dzTerm)
     404       18584 :             dn = Eval_d_Interp_Spl3_pbc((/xs1, xs2, xs3/), TabLR)
     405       18584 :             Nn = SQRT(DOT_PRODUCT(dn, dn))
     406       18584 :             Fterm = Eval_Interp_Spl3_pbc((/xs1, xs2, xs3/), TabLR)
     407        4646 :             tmp1 = ABS(Term - Fterm)
     408       18584 :             tmp2 = SQRT(DOT_PRODUCT(dn - (/dxTerm, dyTerm, dzTerm/), dn - (/dxTerm, dyTerm, dzTerm/)))
     409        4646 :             errf = errf + tmp1
     410        4646 :             maxerrorf = MAX(maxerrorf, tmp1)
     411        4646 :             errd = errd + tmp2
     412        4646 :             maxerrord = MAX(maxerrord, tmp2)
     413             :             WRITE (iw, '(T5,F15.10,5X,F15.10,T50,2F12.9,T80,A,T85,F15.10,5X,F15.10,T130,2F12.9)') &
     414        4646 :                Term, Fterm, tmp1, maxerrorf, "*", Na, Nn, tmp2, maxerrord
     415        4646 :             xs1 = xs1 + dr1
     416        4646 :             xs2 = xs2 + dr2
     417        4692 :             xs3 = xs3 + dr3
     418             :          END DO
     419          46 :          WRITE (iw, '(A,T5,A,T50,F12.9,T130,F12.9)') "#", "Averages", errf/REAL(npoints, kind=dp), &
     420          92 :             errd/REAL(npoints, kind=dp)
     421             :       END IF
     422         102 :       CALL timestop(handle)
     423         102 :       CALL cp_print_key_finished_output(iw, logger, param_section, "check_spline")
     424             : 
     425         102 :    END SUBROUTINE check_spline_interp_TabLR
     426             : 
     427             : END MODULE ewald_spline_util
     428             :