LCOV - code coverage report
Current view: top level - src/pw - pw_grid_info.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.8 % 80 79
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 6 6

            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 This module returns additional info on PW grids
      10              : !> \par History
      11              : !>      JGH (09-06-2007) : Created from pw_grids
      12              : !> \author JGH
      13              : ! **************************************************************************************************
      14              : MODULE pw_grid_info
      15              : 
      16              :    USE fft_tools,                       ONLY: FFT_RADIX_NEXT,&
      17              :                                               FFT_RADIX_NEXT_ODD,&
      18              :                                               fft_radix_operations
      19              :    USE kinds,                           ONLY: dp
      20              :    USE mathconstants,                   ONLY: twopi
      21              :    USE pw_grid_types,                   ONLY: pw_grid_type
      22              : #include "../base/base_uses.f90"
      23              : 
      24              :    IMPLICIT NONE
      25              : 
      26              :    PRIVATE
      27              :    PUBLIC :: pw_find_cutoff, pw_grid_init_setup, pw_grid_bounds_from_n, pw_grid_n_for_fft
      28              : 
      29              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_grid_info'
      30              : 
      31              : CONTAINS
      32              : 
      33              : ! **************************************************************************************************
      34              : !> \brief ...
      35              : !> \param hmat ...
      36              : !> \param cutoff ...
      37              : !> \param spherical ...
      38              : !> \param odd ...
      39              : !> \param fft_usage ...
      40              : !> \param ncommensurate ...
      41              : !> \param icommensurate ...
      42              : !> \param ref_grid ...
      43              : !> \param n_orig ...
      44              : !> \return ...
      45              : ! **************************************************************************************************
      46        35774 :    FUNCTION pw_grid_init_setup(hmat, cutoff, spherical, odd, fft_usage, ncommensurate, &
      47        35774 :                                icommensurate, ref_grid, n_orig) RESULT(n)
      48              : 
      49              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat
      50              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
      51              :       LOGICAL, INTENT(IN)                                :: spherical, odd, fft_usage
      52              :       INTEGER, INTENT(IN)                                :: ncommensurate, icommensurate
      53              :       TYPE(pw_grid_type), INTENT(IN), OPTIONAL           :: ref_grid
      54              :       INTEGER, INTENT(IN), OPTIONAL                      :: n_orig(3)
      55              :       INTEGER, DIMENSION(3)                              :: n
      56              : 
      57              :       INTEGER                                            :: my_icommensurate
      58              : 
      59        35774 :       IF (ncommensurate > 0) THEN
      60         1114 :          my_icommensurate = icommensurate
      61         1114 :          CPASSERT(icommensurate > 0)
      62         1114 :          CPASSERT(icommensurate <= ncommensurate)
      63              :       ELSE
      64              :          my_icommensurate = 0
      65              :       END IF
      66              : 
      67         1114 :       IF (my_icommensurate > 1) THEN
      68          840 :          CPASSERT(PRESENT(ref_grid))
      69         4200 :          n = ref_grid%npts/2**(my_icommensurate - 1)
      70         4200 :          CPASSERT(ALL(ref_grid%npts == n*2**(my_icommensurate - 1)))
      71         3360 :          CPASSERT(ALL(pw_grid_n_for_fft(n) == n))
      72              :       ELSE
      73              :          n = pw_grid_find_n(hmat, cutoff=cutoff, fft_usage=fft_usage, ncommensurate=ncommensurate, &
      74        34934 :                             spherical=spherical, odd=odd, n_orig=n_orig)
      75              :       END IF
      76              : 
      77        35774 :    END FUNCTION pw_grid_init_setup
      78              : 
      79              : ! **************************************************************************************************
      80              : !> \brief returns the n needed for the grid with all the given constraints
      81              : !> \param hmat ...
      82              : !> \param cutoff ...
      83              : !> \param fft_usage ...
      84              : !> \param spherical ...
      85              : !> \param odd ...
      86              : !> \param ncommensurate ...
      87              : !> \param n_orig ...
      88              : !> \return ...
      89              : !> \author fawzi
      90              : ! **************************************************************************************************
      91        34934 :    FUNCTION pw_grid_find_n(hmat, cutoff, fft_usage, spherical, odd, ncommensurate, &
      92        34934 :                            n_orig) RESULT(n)
      93              : 
      94              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat
      95              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
      96              :       LOGICAL, INTENT(IN)                                :: fft_usage, spherical, odd
      97              :       INTEGER, INTENT(IN)                                :: ncommensurate
      98              :       INTEGER, INTENT(IN), OPTIONAL                      :: n_orig(3)
      99              :       INTEGER, DIMENSION(3)                              :: n
     100              : 
     101              :       INTEGER                                            :: idir, my_icommensurate, &
     102              :                                                             my_ncommensurate, nsubgrid, &
     103              :                                                             nsubgrid_new, ntest(3), t_icommensurate
     104              :       LOGICAL                                            :: ftest, subgrid_is_OK
     105              : 
     106              : ! ncommensurate is the number of commensurate grids
     107              : ! in order to have non-commensurate grids ncommensurate must be 0
     108              : ! icommensurte  is the level number of communensurate grids
     109              : ! this implies that the number of grid points in each direction
     110              : ! is k*2**(ncommensurate-icommensurate)
     111              : 
     112        34934 :       my_ncommensurate = ncommensurate
     113        34934 :       IF (my_ncommensurate > 0) THEN
     114              :          my_icommensurate = 1
     115              :       ELSE
     116              :          my_icommensurate = 0
     117              :       END IF
     118        34660 :       CPASSERT(my_icommensurate <= my_ncommensurate)
     119        34934 :       CPASSERT(my_icommensurate > 0 .OR. my_ncommensurate <= 0)
     120        34934 :       CPASSERT(my_ncommensurate >= 0)
     121              : 
     122        34934 :       IF (PRESENT(n_orig)) THEN
     123         9592 :          n = n_orig
     124              :       ELSE
     125        32536 :          CPASSERT(cutoff > 0.0_dp)
     126        32536 :          n = pw_grid_n_from_cutoff(hmat, cutoff)
     127              :       END IF
     128              : 
     129        34934 :       IF (fft_usage) THEN
     130       139736 :          n = pw_grid_n_for_fft(n, odd=odd)
     131              : 
     132        34934 :          IF (.NOT. spherical) THEN
     133       130416 :             ntest = n
     134              : 
     135        32604 :             IF (my_ncommensurate > 0) THEN
     136         1096 :                DO idir = 1, 3
     137          274 :                   DO
     138              :                      ! find valid radix >= ntest
     139         2260 :                      CALL fft_radix_operations(ntest(idir), n(idir), FFT_RADIX_NEXT)
     140              :                      ! check every subgrid of n
     141         2260 :                      subgrid_is_OK = .TRUE.
     142         4780 :                      DO t_icommensurate = 1, my_ncommensurate - 1
     143         3958 :                         nsubgrid = n(idir)/2**(my_ncommensurate - t_icommensurate)
     144         3958 :                         CALL fft_radix_operations(nsubgrid, nsubgrid_new, FFT_RADIX_NEXT)
     145              :                         subgrid_is_OK = (nsubgrid == nsubgrid_new) .AND. &
     146         3958 :                                         (MODULO(n(idir), 2**(my_ncommensurate - t_icommensurate)) == 0)
     147          822 :                         IF (.NOT. subgrid_is_OK) EXIT
     148              :                      END DO
     149         2260 :                      IF (subgrid_is_OK) THEN
     150              :                         EXIT
     151              :                      ELSE
     152              :                         ! subgrid wasn't OK, increment ntest and try again
     153         1438 :                         ntest(idir) = n(idir) + 1
     154              :                      END IF
     155              :                   END DO
     156              :                END DO
     157              :             END IF
     158              :          END IF
     159              :       ELSE
     160              :          ! without a cutoff and HALFSPACE we have to be sure that there is
     161              :          ! a negative counterpart to every g vector (-> odd number of grid points)
     162            0 :          IF (odd) n = n + MOD(n + 1, 2)
     163              : 
     164              :       END IF
     165              : 
     166              :       ! final check if all went fine ...
     167         2604 :       IF (my_ncommensurate > 0) THEN
     168         1388 :          DO my_icommensurate = 1, my_ncommensurate
     169         5570 :             ftest = ANY(MODULO(n, 2**(my_ncommensurate - my_icommensurate)) .NE. 0)
     170         1388 :             CPASSERT(.NOT. ftest)
     171              :          END DO
     172              :       END IF
     173              : 
     174        34934 :    END FUNCTION pw_grid_find_n
     175              : 
     176              : ! **************************************************************************************************
     177              : !> \brief returns the closest number of points >= n, on which you can perform
     178              : !>      ffts
     179              : !> \param n the minimum number of points you want
     180              : !> \param odd if the number has to be odd
     181              : !> \return ...
     182              : !> \author fawzi
     183              : !> \note
     184              : !>      result<=n
     185              : ! **************************************************************************************************
     186        36210 :    FUNCTION pw_grid_n_for_fft(n, odd) RESULT(nout)
     187              :       INTEGER, DIMENSION(3), INTENT(in)                  :: n
     188              :       LOGICAL, INTENT(in), OPTIONAL                      :: odd
     189              :       INTEGER, DIMENSION(3)                              :: nout
     190              : 
     191              :       LOGICAL                                            :: my_odd
     192              : 
     193        36210 :       my_odd = .FALSE.
     194        36210 :       IF (PRESENT(odd)) my_odd = odd
     195       144840 :       CPASSERT(ALL(n >= 0))
     196        36210 :       IF (my_odd) THEN
     197          468 :          CALL fft_radix_operations(n(1), nout(1), FFT_RADIX_NEXT_ODD)
     198          468 :          CALL fft_radix_operations(n(2), nout(2), FFT_RADIX_NEXT_ODD)
     199          468 :          CALL fft_radix_operations(n(3), nout(3), FFT_RADIX_NEXT_ODD)
     200              :       ELSE
     201        35742 :          CALL fft_radix_operations(n(1), nout(1), FFT_RADIX_NEXT)
     202        35742 :          CALL fft_radix_operations(n(2), nout(2), FFT_RADIX_NEXT)
     203        35742 :          CALL fft_radix_operations(n(3), nout(3), FFT_RADIX_NEXT)
     204              :       END IF
     205              : 
     206        36210 :    END FUNCTION pw_grid_n_for_fft
     207              : 
     208              : ! **************************************************************************************************
     209              : !> \brief Find the number of points that give at least the requested cutoff
     210              : !> \param hmat ...
     211              : !> \param cutoff ...
     212              : !> \return ...
     213              : !> \par History
     214              : !>      JGH (21-12-2000) : Simplify parameter list, bounds will be global
     215              : !>      JGH ( 8-01-2001) : Add check to FFT allowd grids (this now depends
     216              : !>                         on the FFT library.
     217              : !>                         Should the pw_grid_type have a reference to the FFT
     218              : !>                         library ?
     219              : !>      JGH (28-02-2001) : Only do conditional check for FFT
     220              : !>      JGH (21-05-2002) : Optimise code, remove orthorhombic special case
     221              : !> \author apsi
     222              : !>      Christopher Mundy
     223              : ! **************************************************************************************************
     224        32536 :    FUNCTION pw_grid_n_from_cutoff(hmat, cutoff) RESULT(n)
     225              : 
     226              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: hmat
     227              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
     228              :       INTEGER, DIMENSION(3)                              :: n
     229              : 
     230              :       INTEGER                                            :: i
     231              :       REAL(KIND=dp)                                      :: alat(3)
     232              : 
     233       130144 :       DO i = 1, 3
     234       422968 :          alat(i) = SUM(hmat(:, i)**2)
     235              :       END DO
     236       130144 :       CPASSERT(ALL(alat /= 0._dp))
     237       130144 :       n = 2*FLOOR(SQRT(2.0_dp*cutoff*alat)/twopi) + 1
     238              : 
     239        32536 :    END FUNCTION pw_grid_n_from_cutoff
     240              : 
     241              : ! **************************************************************************************************
     242              : !> \brief returns the bounds that distribute n points evenly around 0
     243              : !> \param npts the number of points in each direction
     244              : !> \return ...
     245              : !> \author fawzi
     246              : ! **************************************************************************************************
     247        33426 :    FUNCTION pw_grid_bounds_from_n(npts) RESULT(bounds)
     248              :       INTEGER, DIMENSION(3), INTENT(in)                  :: npts
     249              :       INTEGER, DIMENSION(2, 3)                           :: bounds
     250              : 
     251       133704 :       bounds(1, :) = -npts/2
     252       133704 :       bounds(2, :) = bounds(1, :) + npts - 1
     253              : 
     254        33426 :    END FUNCTION pw_grid_bounds_from_n
     255              : 
     256              : ! **************************************************************************************************
     257              : !> \brief Given a grid and a box, calculate the corresponding cutoff
     258              : !>      *** This routine calculates the cutoff in MOMENTUM UNITS! ***
     259              : !> \param npts ...
     260              : !> \param h_inv ...
     261              : !> \return ...
     262              : !> \par History
     263              : !>      JGH (20-12-2000) : Deleted some strange comments
     264              : !> \author apsi
     265              : !>      Christopher Mundy
     266              : !> \note
     267              : !>      This routine is local. It works independent from the distribution
     268              : !>      of PW on processors.
     269              : !>      npts is the grid size for the full box.
     270              : ! **************************************************************************************************
     271         3409 :    FUNCTION pw_find_cutoff(npts, h_inv) RESULT(cutoff)
     272              : 
     273              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: npts
     274              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: h_inv
     275              :       REAL(KIND=dp)                                      :: cutoff
     276              : 
     277              :       REAL(KIND=dp)                                      :: gcut, gdum(3), length
     278              : 
     279              : ! compute 2*pi*h_inv^t*g  where g = (nmax[1],0,0)
     280              : 
     281        13636 :       gdum(:) = twopi*h_inv(1, :)*REAL((npts(1) - 1)/2, KIND=dp)
     282         3409 :       length = SQRT(gdum(1)**2 + gdum(2)**2 + gdum(3)**2)
     283         3409 :       gcut = length
     284              : 
     285              :       ! compute 2*pi*h_inv^t*g  where g = (0,nmax[2],0)
     286        13636 :       gdum(:) = twopi*h_inv(2, :)*REAL((npts(2) - 1)/2, KIND=dp)
     287         3409 :       length = SQRT(gdum(1)**2 + gdum(2)**2 + gdum(3)**2)
     288         3409 :       gcut = MIN(gcut, length)
     289              : 
     290              :       ! compute 2*pi*h_inv^t*g  where g = (0,0,nmax[3])
     291        13636 :       gdum(:) = twopi*h_inv(3, :)*REAL((npts(3) - 1)/2, KIND=dp)
     292         3409 :       length = SQRT(gdum(1)**2 + gdum(2)**2 + gdum(3)**2)
     293         3409 :       gcut = MIN(gcut, length)
     294              : 
     295         3409 :       cutoff = gcut - 1.e-8_dp
     296              : 
     297         3409 :    END FUNCTION pw_find_cutoff
     298              : 
     299              : END MODULE pw_grid_info
     300              : 
        

Generated by: LCOV version 2.0-1