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