pw_unittest.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
PROGRAM pw_unittest
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE kinds,                           ONLY: dp
   USE libcp2k,                         ONLY: cp2k_finalize,&
                                              cp2k_init
   USE mathconstants,                   ONLY: fourpi,&
                                              pi
   USE message_passing,                 ONLY: mp_comm_type,&
                                              mp_comm_world,&
                                              mp_para_env_type
   USE mt_util,                         ONLY: mt2d_kernel_strain
   USE pw_copy_all,                     ONLY: pw_copy_match
   USE pw_grid_types,                   ONLY: fullspace,&
                                              halfspace,&
                                              pw_grid_type
   USE pw_grids,                        ONLY: pw_grid_create,&
                                              pw_grid_release
   USE pw_methods,                      ONLY: do_accurate_sum,&
                                              do_standard_sum,&
                                              pw_integral_ab
   USE pw_poisson_methods,              ONLY: pw_poisson_set,&
                                              pw_poisson_solve
   USE pw_poisson_types,                ONLY: pw_poisson_mt,&
                                              pw_poisson_type
   USE pw_pool_types,                   ONLY: pw_pool_create
   USE pw_types,                        ONLY: pw_c1d_gs_type,&
                                              pw_r1d_gs_type
   USE realspace_grid_cube,             ONLY: cube_read_values,&
                                              cube_value_format,&
                                              cube_values_format
 
   IMPLICIT NONE
 
   INTEGER                                            :: mode
   REAL(dp)                                           :: max_error
   TYPE(cp_logger_type), POINTER                      :: logger
   TYPE(mp_para_env_type), POINTER                     :: para_env
 
   CALL cp2k_init()
   logger => cp_get_default_logger()
   para_env => logger%para_env
   CALL check_cube_values()
   DO mode = 0, 8
      CALL check_copy(mode)
   END DO
   max_error = 0.0_dp
   CALL check_pw_integral_and_poisson()
   CALL check_mt2d_kernel()
   print '(A,ES16.7)', 'PW integral and MT2D kernel maximum scaled error: ', max_error
   CALL cp2k_finalize()
 
CONTAINS
 
! **************************************************************************************************
!> \brief Preserve three-digit cube exponents and parse adjacent values across z-slices.
! **************************************************************************************************
   SUBROUTINE check_cube_values()
      CHARACTER(LEN=13)                                  :: value
      CHARACTER(LEN=79)                                  :: values
      CHARACTER(LEN=92)                                  :: slice_values
      REAL(KIND=dp), DIMENSION(6)                        :: buffer, reference
      REAL(KIND=dp), DIMENSION(7)                        :: slice_buffer, slice_reference
 
      WRITE (value, cube_value_format) 3.33004e-101_dp
      IF (index(value, "E-101") <= 0) THEN
         error stop "Cube value format must preserve three-digit negative exponents."
      END IF
 
      reference = [0.0_dp, -0.91458e-35_dp, 3.33004e-101_dp, &
                   -1.32487e-103_dp, 1.0_dp, -1.0_dp]
      WRITE (values(1:78), cube_values_format) reference
      values(79:79) = new_line("C")
      IF (index(values, "E-101") <= 0) THEN
         error stop "Cube line format must preserve E-101."
      END IF
      IF (index(values, "E-103") <= 0) THEN
         error stop "Cube line format must preserve E-103."
      END IF
      CALL cube_read_values(values, buffer)
      IF (maxval(abs(buffer - reference)) > 1.0e-12_dp) THEN
         error stop "Cube reader must parse adjacent explicit-exponent values."
      END IF
 
      slice_reference(1:6) = reference
      slice_reference(7) = 0.27183e-123_dp
      WRITE (slice_values(1:78), cube_values_format) slice_reference(1:6)
      slice_values(79:79) = new_line("C")
      WRITE (slice_values(80:92), cube_value_format) slice_reference(7)
      CALL cube_read_values(slice_values, slice_buffer)
      IF (maxval(abs(slice_buffer - slice_reference)) > 1.0e-12_dp) THEN
         error stop "Cube reader must parse multi-line z-slices."
      END IF
   END SUBROUTINE check_cube_values
 
! **************************************************************************************************
!> \brief Compare identical, permuted, redistributed, partial and empty PW grids.
!> \param mode Layout variant; mode two changes only one rank's ordering.
!> \author Thomas D. Kuehne (tkuehne@cp2k.org)
! **************************************************************************************************
   SUBROUTINE check_copy(mode)
      INTEGER, INTENT(IN)                                :: mode
 
      COMPLEX(dp), PARAMETER :: untouched = (-999.0_dp, 17.0_dp)
 
      COMPLEX(dp)                                        :: expected
      INTEGER                                            :: dimensions(1), i, j, nsource, ntarget, &
                                                            owner, rank
      TYPE(pw_c1d_gs_type)                               :: destination, source
      TYPE(pw_grid_type), TARGET                         :: source_grid, target_grid
 
      rank = para_env%mepos
      nsource = 7
      IF (mode == 6 .OR. (mode >= 7 .AND. rank == 0)) nsource = 0
      owner = rank
      IF (mode == 3 .OR. mode == 8) owner = mod(rank + 1, para_env%num_pe)
      ntarget = nsource
      IF (mode == 4) ntarget = nsource - 1
      IF (mode == 5) ntarget = nsource + 1
      IF (mode == 8) ntarget = merge(0, 7, owner == 0)
      ALLOCATE (source_grid%g_hat(3, nsource), target_grid%g_hat(3, ntarget))
      ALLOCATE (source%array(nsource), destination%array(ntarget))
      source%pw_grid => source_grid
      destination%pw_grid => target_grid
      dimensions = [para_env%num_pe]
      CALL source_grid%para%group%create(para_env, 1, dimensions)
      target_grid%para%group = source_grid%para%group
      source_grid%ngpts_cut_local = nsource
      target_grid%ngpts_cut_local = ntarget
      DO i = 1, nsource
         source_grid%g_hat(:, i) = [i, rank, 0]
         source%array(i) = cmplx(i + 11*rank, -i - 3*rank, dp)
      END DO
      DO i = 1, ntarget
         j = i
         IF (mode == 1 .OR. (mode == 2 .AND. rank == 0)) j = ntarget + 1 - i
         target_grid%g_hat(:, i) = [j, owner, 0]
      END DO
      destination%array(:) = untouched
      CALL pw_copy_match(source, destination)
      DO i = 1, ntarget
         j = target_grid%g_hat(1, i)
         expected = cmplx(j + 11*owner, -j - 3*owner, dp)
         IF (mode == 5 .AND. j > nsource) expected = untouched
         IF (abs(destination%array(i) - expected) > 1.e-14_dp) error stop 'Incorrect PW copy'
      END DO
      DEALLOCATE (source%array, destination%array)
      CALL source_grid%para%group%free()
   END SUBROUTINE check_copy
 
! **************************************************************************************************
!> \brief Check PW multiplicity and volume normalization, and auxiliary-density Poisson routes.
! **************************************************************************************************
   SUBROUTINE check_pw_integral_and_poisson()
      REAL(dp), PARAMETER :: step = 5.e-4_dp, &
         strain_steps(4) = [step, -step, 2.0_dp*step, -2.0_dp*step]
 
      COMPLEX(dp)                                        :: a_mode, b_mode
      INTEGER                                            :: axis, ig, sign_step, span, sumtype, test
      REAL(dp) :: cell(3, 3), derivative, energy, energy_dv, energy_v, energy_v_dv, expected, &
         expected_real, kernel, reference, shifted_g(3), strain_energy(4), stress(3, 3), &
         stress_dv(3, 3), stress_only(3, 3), stress_v(3, 3), stress_v_dv(3, 3)
      TYPE(mp_comm_type)                                 :: comm
      TYPE(pw_c1d_gs_type)                               :: a, b, gradient(3), gradient_v_dv(3), &
                                                            potential, potential_v_dv
      TYPE(pw_grid_type), POINTER                        :: grid
      TYPE(pw_poisson_type)                              :: poisson
      TYPE(pw_r1d_gs_type)                               :: a_real, b_real
 
      CALL comm%set_handle(mp_comm_world%get_handle())
      ! Nonunit volume and an appreciable truncation factor exp(-pi).
      cell = 0.0_dp
      DO axis = 1, 3
         cell(axis, axis) = 7.0_dp**(1.0_dp/3.0_dp)
      END DO
      a_mode = cmplx(0.4_dp, -0.2_dp, dp)
      b_mode = cmplx(-0.7_dp, 0.9_dp, dp)
      NULLIFY (grid)
      DO span = fullspace, halfspace, -1
         ! Reciprocal inputs and outputs need no FFT. The production grid still
         ! supplies the distributed G0 ownership and half/full-space mappings.
         CALL pw_grid_create(grid, comm, cell, grid_span=span, npts=[5, 5, 5], odd=.true., fft_usage=.false.)
         CALL a%create(grid)
         CALL b%create(grid)
         CALL a_real%create(grid)
         CALL b_real%create(grid)
         CALL potential%create(grid)
         CALL potential_v_dv%create(grid)
         DO axis = 1, 3
            CALL gradient(axis)%create(grid)
            CALL gradient_v_dv(axis)%create(grid)
         END DO
         ALLOCATE (poisson%pw_pools(1))
         CALL pw_pool_create(poisson%pw_pools(1)%pool, grid)
         poisson%pw_level = 1
         poisson%parameters%solver = pw_poisson_mt
         poisson%parameters%periodic = [1, 1, 0]
         poisson%parameters%mt_alpha = 1.0_dp
         CALL pw_poisson_set(poisson, cell_hmat=cell, use_level=1)
         DO test = 1, 2
            ! Constant G0 only and G0 plus a Hermitian mode pair.
            a%array = 0.0_dp
            b%array = 0.0_dp
            expected = 0.0_dp
            expected_real = 0.0_dp
            reference = 0.0_dp
            IF (grid%have_g0) THEN
               a%array(1) = 2.0_dp
               b%array(1) = 3.0_dp
            END IF
            expected = 6.0_dp
            expected_real = 6.0_dp
            reference = 6.0_dp*mt2d_kernel_reference([0.0_dp, 0.0_dp, 0.0_dp], 3, cell(3, 3))
            IF (test /= 1) THEN
               DO ig = 1, grid%ngpts_cut_local
                  IF (any(grid%g_hat(2:3, ig) /= 0)) cycle
                  IF (grid%g_hat(1, ig) == 1) THEN
                     a%array(ig) = a_mode
                     b%array(ig) = b_mode
                  ELSE IF (grid%g_hat(1, ig) == -1) THEN
                     a%array(ig) = conjg(a_mode)
                     b%array(ig) = conjg(b_mode)
                  END IF
               END DO
               expected = expected + 2.0_dp*real(conjg(a_mode)*b_mode, dp)
               expected_real = expected_real + 2.0_dp*real(a_mode, dp)*real(b_mode, dp)
               kernel = mt2d_kernel_reference([2.0_dp*pi/cell(1, 1), 0.0_dp, 0.0_dp], 3, cell(3, 3))
               reference = reference + 2.0_dp*real(conjg(a_mode)*b_mode, dp)*kernel
            END IF
            a_real%array = real(a%array, dp)
            b_real%array = real(b%array, dp)
            ! G0 must have weight one, the mode pair weight two, and both must
            ! receive V=7 for an integral. just_sum retains the unscaled oracle.
            DO sumtype = do_accurate_sum, do_standard_sum
               CALL check_error(pw_integral_ab(a, b, sumtype=sumtype), 7.0_dp*expected, 3.e-14_dp)
               CALL check_error(pw_integral_ab(a, b, sumtype=sumtype, just_sum=.true.), expected, 3.e-14_dp)
               CALL check_error(pw_integral_ab(a_real, b_real, sumtype=sumtype), 7.0_dp*expected_real, 3.e-14_dp)
               CALL check_error(pw_integral_ab(a_real, b_real, sumtype=sumtype, just_sum=.true.), &
                                expected_real, 3.e-14_dp)
               CALL check_error(pw_integral_ab(a, b_real, sumtype=sumtype), 7.0_dp*expected_real, 3.e-14_dp)
               CALL check_error(pw_integral_ab(a_real, b, sumtype=sumtype), 7.0_dp*expected_real, 3.e-14_dp)
            END DO
            CALL pw_poisson_solve(poisson, a, energy, h_stress=stress, aux_density=b)
            CALL pw_poisson_solve(poisson, a, energy_dv, dvhartree=gradient, h_stress=stress_dv, aux_density=b)
            CALL pw_poisson_solve(poisson, a, energy_v, potential, h_stress=stress_v, aux_density=b)
            CALL pw_poisson_solve(poisson, a, energy_v_dv, potential_v_dv, dvhartree=gradient_v_dv, &
                                  h_stress=stress_v_dv, aux_density=b)
            CALL pw_poisson_solve(poisson, a, h_stress=stress_only, aux_density=b)
            CALL check_error(energy, 0.5_dp*7.0_dp*reference, 3.e-13_dp)
            CALL check_error(energy_dv, energy, 3.e-13_dp)
            CALL check_error(energy_v, energy, 3.e-13_dp)
            CALL check_error(energy_v_dv, energy, 3.e-13_dp)
            CALL check_error(maxval(abs(stress_dv - stress)), 0.0_dp, 3.e-11_dp)
            CALL check_error(maxval(abs(stress_v - stress)), 0.0_dp, 3.e-11_dp)
            CALL check_error(maxval(abs(stress_v_dv - stress)), 0.0_dp, 3.e-11_dp)
            CALL check_error(maxval(abs(stress_only - stress)), 0.0_dp, 3.e-11_dp)
            IF (test == 2) THEN
               ! Fixed Fourier charges Q=V*rho give E(epsilon)=V/2/(1+epsilon)
               ! times the strained kernel sum. This absolute production-virial
               ! oracle includes both G0/volume and the nonzero-mode kernel term.
               DO sign_step = 1, 4
                  shifted_g = [2.0_dp*pi/cell(1, 1)/(1.0_dp + strain_steps(sign_step)), 0.0_dp, 0.0_dp]
                  reference = 6.0_dp*mt2d_kernel_reference([0.0_dp, 0.0_dp, 0.0_dp], 3, cell(3, 3)) &
                              + 2.0_dp*real(conjg(a_mode)*b_mode, dp)*mt2d_kernel_reference(shifted_g, 3, cell(3, 3))
                  strain_energy(sign_step) = 0.5_dp*7.0_dp*reference/(1.0_dp + strain_steps(sign_step))
               END DO
               derivative = (8.0_dp*(strain_energy(1) - strain_energy(2)) &
                             - strain_energy(3) + strain_energy(4))/(12.0_dp*step)
               CALL check_error(stress(1, 1) + derivative, 0.0_dp, 3.e-9_dp)
            END IF
            DO ig = 1, grid%ngpts_cut_local
               kernel = mt2d_kernel_reference(grid%g(:, ig), 3, cell(3, 3))
               CALL check_error(abs(potential%array(ig) - kernel*a%array(ig)), 0.0_dp, 3.e-13_dp)
               CALL check_error(abs(potential_v_dv%array(ig) - kernel*a%array(ig)), 0.0_dp, 3.e-13_dp)
               DO axis = 1, 3
                  CALL check_error(abs(gradient(axis)%array(ig) &
                                       - cmplx(0.0_dp, grid%g(axis, ig), dp)*kernel*a%array(ig)), 0.0_dp, 3.e-13_dp)
                  CALL check_error(abs(gradient_v_dv(axis)%array(ig) &
                                       - cmplx(0.0_dp, grid%g(axis, ig), dp)*kernel*a%array(ig)), 0.0_dp, 3.e-13_dp)
               END DO
            END DO
         END DO
         CALL poisson%release()
         CALL a%release()
         CALL b%release()
         CALL a_real%release()
         CALL b_real%release()
         CALL potential%release()
         CALL potential_v_dv%release()
         DO axis = 1, 3
            CALL gradient(axis)%release()
            CALL gradient_v_dv(axis)%release()
         END DO
         CALL pw_grid_release(grid)
      END DO
   END SUBROUTINE check_pw_integral_and_poisson
 
! **************************************************************************************************
!> \brief Compare the MT2D kernel and in-plane virial with independent strained-cell energies.
! **************************************************************************************************
   SUBROUTINE check_mt2d_kernel()
      REAL(dp), PARAMETER :: step = 3.e-4_dp, shifts(4) = [step, -step, 2.0_dp*step, -2.0_dp*step]
 
      INTEGER                                            :: alpha, beta, ip, iz, jp, miller(3), &
                                                            mode, modes(3, 5), sign_step
      REAL(dp)                                           :: derivative, g(3), h_inv(3, 3), kernel, &
                                                            lengths(3), shifted_g(3), &
                                                            shifted_kernel(4), strain(3, 3)
 
      ! Rotate the Cartesian vacuum axis. A shear within the periodic plane is
      ! permitted; only normal/periodic mixing is excluded by the MT2D kernel.
      DO iz = 1, 3
         ip = mod(iz, 3) + 1
         jp = mod(iz + 1, 3) + 1
         lengths = [9.0_dp, 13.0_dp, 17.0_dp]
         h_inv = 0.0_dp
         DO alpha = 1, 3
            h_inv(alpha, alpha) = 1.0_dp/lengths(alpha)
         END DO
         h_inv(ip, jp) = -1.4_dp/(lengths(ip)*lengths(jp))
         modes = 0
         ! G=0, both odd and even normal-only modes, and several modes with
         ! appreciable exponential truncation reject the periodic K**2 formula.
         modes(iz, 2:3) = [1, 2]
         modes(ip, 4) = 1
         modes(:, 5) = [1, 2, 1]
         DO mode = 1, SIZE(modes, 2)
            miller = modes(:, mode)
            g = 2.0_dp*pi*matmul(transpose(h_inv), real(miller, dp))
            CALL mt2d_kernel_strain(g, iz, lengths(iz), kernel, strain)
            CALL check_error(kernel, mt2d_kernel_reference(g, iz, lengths(iz)), 3.e-14_dp)
            CALL check_error(maxval(abs(strain(iz, :))), 0.0_dp, 0.0_dp)
            CALL check_error(maxval(abs(strain(:, iz))), 0.0_dp, 0.0_dp)
            IF (mode <= 3) CALL check_error(maxval(abs(strain)), 0.0_dp, 0.0_dp)
            DO alpha = 1, 3
               IF (alpha == iz) cycle
               DO beta = 1, 3
                  IF (beta == iz) cycle
                  DO sign_step = 1, 4
                     ! g'=(I+epsilon_ab)^(-T) g, evaluated independently of dK.
                     shifted_g = g
                     IF (alpha == beta) THEN
                        shifted_g(beta) = g(beta)/(1.0_dp + shifts(sign_step))
                     ELSE
                        shifted_g(beta) = g(beta) - shifts(sign_step)*g(alpha)
                     END IF
                     shifted_kernel(sign_step) = mt2d_kernel_reference(shifted_g, iz, lengths(iz))
                  END DO
                  derivative = (8.0_dp*(shifted_kernel(1) - shifted_kernel(2)) &
                                - shifted_kernel(3) + shifted_kernel(4))/(12.0_dp*step)
                  CALL check_error(strain(alpha, beta), derivative, 3.e-9_dp)
               END DO
            END DO
         END DO
      END DO
   END SUBROUTINE check_mt2d_kernel
 
! **************************************************************************************************
!> \brief Evaluate the original MT2D influence kernel independently of its strain derivative.
!> \param g ...
!> \param iz ...
!> \param slab_size ...
!> \return ...
! **************************************************************************************************
   FUNCTION mt2d_kernel_reference(g, iz, slab_size) RESULT(kernel)
      REAL(dp), INTENT(IN)                               :: g(3)
      INTEGER, INTENT(IN)                                :: iz
      REAL(dp), INTENT(IN)                               :: slab_size
      REAL(dp)                                           :: kernel
 
      REAL(dp)                                           :: g2, g_parallel(3), radius
 
      g2 = sum(g**2)
      IF (g2 == 0.0_dp) THEN
         kernel = pi*slab_size**2/2.0_dp
      ELSE
         g_parallel = g
         g_parallel(iz) = 0.0_dp
         radius = slab_size/2.0_dp
         kernel = fourpi/g2 - fourpi/g2*cos(g(iz)*radius)*exp(-norm2(g_parallel)*radius)
      END IF
   END FUNCTION mt2d_kernel_reference
 
! **************************************************************************************************
!> \brief Require the scaled error to be finite and within the supplied tolerance.
!> \param value ...
!> \param reference ...
!> \param tolerance ...
! **************************************************************************************************
   SUBROUTINE check_error(value, reference, tolerance)
      REAL(dp), INTENT(IN)                               :: value, reference, tolerance
 
      REAL(dp)                                           :: error
 
      error = abs(value - reference)/max(1.0_dp, abs(reference))
      max_error = max(max_error, error)
      IF (.NOT. (error <= tolerance)) THEN
         print *, 'Coulomb check: value, reference, error, tolerance:', value, reference, error, tolerance
         error stop 'PW integral or MT2D kernel check failed'
      END IF
   END SUBROUTINE check_error
 
END PROGRAM pw_unittest