Line data    Source code

       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8            2 : PROGRAM cell_types_TEST
       9              :    USE cell_types,                      ONLY: cell_type,&
      10              :                                               pbc_stable,&
      11            2 :                                               real_to_scaled,&
      12              :                                               scaled_to_real
      13              :    USE kinds,                           ONLY: dp
      14              : 
      15              :    IMPLICIT NONE
      16              : 
      17              :    TYPE(cell_type), POINTER                           :: cell
      18              : 
      19          138 :    ALLOCATE (cell)
      20           26 :    cell%hmat = 0.0_dp
      21           26 :    cell%h_inv = 0.0_dp
      22            2 :    cell%hmat(1, 1) = 2.0_dp
      23            2 :    cell%hmat(2, 2) = 3.0_dp
      24            2 :    cell%hmat(3, 3) = 4.0_dp
      25            2 :    cell%h_inv(1, 1) = 0.5_dp
      26            2 :    cell%h_inv(2, 2) = 1.0_dp/3.0_dp
      27            2 :    cell%h_inv(3, 3) = 0.25_dp
      28            8 :    cell%perd = 1
      29              : 
      30            2 :    cell%orthorhombic = .TRUE.
      31            2 :    CALL check_half_cell_boundaries(cell)
      32              : 
      33            2 :    cell%orthorhombic = .FALSE.
      34            2 :    cell%hmat(1, 2) = 0.5_dp
      35            2 :    cell%h_inv(1, 2) = -1.0_dp/12.0_dp
      36            2 :    CALL check_half_cell_boundaries(cell)
      37              : 
      38            2 :    DEALLOCATE (cell)
      39              : 
      40              : CONTAINS
      41              : 
      42              : ! **************************************************************************************************
      43              : !> \brief Check that roundoff around half-cell boundaries does not change the periodic image.
      44              : !> \param cell ...
      45              : ! **************************************************************************************************
      46            4 :    SUBROUTINE check_half_cell_boundaries(cell)
      47              : 
      48              :       TYPE(cell_type), POINTER                           :: cell
      49              : 
      50              :       INTEGER, DIMENSION(3), PARAMETER :: lattice_translation = [-1, 1, -2]
      51              :       REAL(KIND=dp), PARAMETER                           :: delta = 8.0_dp*EPSILON(1.0_dp)
      52              : 
      53              :       REAL(KIND=dp), DIMENSION(3)                        :: r, scaled, wrapped, wrapped_scaled, &
      54              :                                                             wrapped_translated, &
      55              :                                                             wrapped_translated_scaled
      56              : 
      57            4 :       scaled = [0.5_dp - delta, -0.5_dp + delta, 1.5_dp - delta]
      58            4 :       CALL scaled_to_real(r, scaled, cell)
      59            4 :       wrapped = pbc_stable(r, cell)
      60            4 :       CALL real_to_scaled(wrapped_scaled, wrapped, cell)
      61           16 :       IF (ANY(SIGN(1.0_dp, wrapped_scaled) /= [-1.0_dp, -1.0_dp, -1.0_dp])) THEN
      62            0 :          ERROR STOP "Roundoff changed the image at a half-cell boundary"
      63              :       END IF
      64           16 :       IF (ANY(ABS(ABS(wrapped_scaled) - 0.5_dp) > 2.0_dp*delta)) THEN
      65            0 :          ERROR STOP "Unexpected wrapped coordinate at a half-cell boundary"
      66              :       END IF
      67              : 
      68           16 :       CALL scaled_to_real(r, scaled + REAL(lattice_translation, KIND=dp), cell)
      69            4 :       wrapped_translated = pbc_stable(r, cell)
      70            4 :       CALL real_to_scaled(wrapped_translated_scaled, wrapped_translated, cell)
      71           16 :       IF (ANY(ABS(wrapped_translated_scaled - wrapped_scaled) > 4.0_dp*delta)) THEN
      72            0 :          ERROR STOP "Periodic wrapping changed under a lattice translation"
      73              :       END IF
      74              : 
      75              :       scaled = [0.5_dp - 128.0_dp*EPSILON(1.0_dp), &
      76              :                 -0.5_dp - 128.0_dp*EPSILON(1.0_dp), &
      77            4 :                 1.5_dp + 128.0_dp*EPSILON(1.0_dp)]
      78            4 :       CALL scaled_to_real(r, scaled, cell)
      79            4 :       wrapped = pbc_stable(r, cell)
      80            4 :       CALL real_to_scaled(wrapped_scaled, wrapped, cell)
      81           16 :       IF (ANY(SIGN(1.0_dp, wrapped_scaled) /= [1.0_dp, 1.0_dp, -1.0_dp])) THEN
      82            0 :          ERROR STOP "Coordinate outside the boundary tolerance changed image"
      83              :       END IF
      84              : 
      85            4 :    END SUBROUTINE check_half_cell_boundaries
      86              : 
      87              : END PROGRAM cell_types_TEST