1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief for a given dr()/dh(r) this will provide the bounds to be used if 8!> one wants to go over a sphere-subregion of given radius 9!> \note 10!> the computation of the exact sphere radius is sensitive to roundoff (e.g. 11!> compiler optimization level) and hence this small roundoff can result in 12!> energy difference of about EPS_DEFAULT in QS energies (one gridpoint more or 13!> less in the density mapping) 14!> \author Joost VandeVondele 15! ************************************************************************************************** 16MODULE cube_utils 17 18 USE kinds, ONLY: dp 19 USE mathlib, ONLY: matvec_3x3 20 USE realspace_grid_types, ONLY: realspace_grid_desc_type 21#include "../base/base_uses.f90" 22 23 IMPLICIT NONE 24 PRIVATE 25 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cube_utils' 26 27 PUBLIC :: cube_info_type 28 29 PUBLIC :: init_cube_info, destroy_cube_info, & 30 return_cube, return_cube_max_iradius, return_cube_nonortho, & 31 compute_cube_center 32 33 TYPE :: cube_ptr 34 INTEGER, POINTER, DIMENSION(:) :: p 35 END TYPE cube_ptr 36 37 TYPE :: cube_info_type 38 INTEGER :: max_radius 39 REAL(KIND=dp) :: dr(3), drmin 40 REAL(KIND=dp) :: dh(3, 3) 41 REAL(KIND=dp) :: dh_inv(3, 3) 42 LOGICAL :: orthorhombic 43 INTEGER, POINTER :: lb_cube(:, :) 44 INTEGER, POINTER :: ub_cube(:, :) 45 TYPE(cube_ptr), POINTER, DIMENSION(:) :: sphere_bounds 46 INTEGER, POINTER :: sphere_bounds_count(:) 47 REAL(KIND=dp) :: max_rad_ga 48 END TYPE cube_info_type 49 50CONTAINS 51! ************************************************************************************************** 52!> \brief unifies the computation of the cube center, so that differences in 53!> implementation, and thus roundoff and numerics can not lead to 54!> off-by-one errors (which can lead to out-of-bounds access with distributed grids). 55!> in principle, something similar would be needed for the computation of the cube bounds 56!> 57!> \param cube_center ... 58!> \param rs_desc ... 59!> \param zeta ... 60!> \param zetb ... 61!> \param ra ... 62!> \param rab ... 63!> \par History 64!> 11.2008 created [Joost VandeVondele] 65! ************************************************************************************************** 66 SUBROUTINE compute_cube_center(cube_center, rs_desc, zeta, zetb, ra, rab) 67 68 INTEGER, DIMENSION(3), INTENT(OUT) :: cube_center 69 TYPE(realspace_grid_desc_type), POINTER :: rs_desc 70 REAL(KIND=dp), INTENT(IN) :: zeta, zetb, ra(3), rab(3) 71 72 REAL(KIND=dp) :: zetp 73 REAL(KIND=dp), DIMENSION(3) :: rp 74 75 zetp = zeta + zetb 76 rp(:) = ra(:) + zetb/zetp*rab(:) 77 cube_center(:) = FLOOR(MATMUL(rs_desc%dh_inv, rp)) 78 79 END SUBROUTINE compute_cube_center 80 81! ************************************************************************************************** 82!> \brief ... 83!> \param info ... 84!> \param radius ... 85!> \param lb ... 86!> \param ub ... 87!> \param rp ... 88! ************************************************************************************************** 89 SUBROUTINE return_cube_nonortho(info, radius, lb, ub, rp) 90 91 TYPE(cube_info_type), INTENT(IN) :: info 92 REAL(KIND=dp), INTENT(IN) :: radius 93 INTEGER, INTENT(OUT) :: lb(3), ub(3) 94 REAL(KIND=dp), INTENT(IN) :: rp(3) 95 96 CHARACTER(LEN=*), PARAMETER :: routineN = 'return_cube_nonortho', & 97 routineP = moduleN//':'//routineN 98 99 INTEGER :: i, j, k 100 REAL(KIND=dp) :: point(3), res(3) 101 102 IF (radius > info%max_rad_ga) THEN 103 ! 104 ! This is an important check. If the required radius for mapping the density is different 105 ! from the actual computed one, (significant) errors can occur. 106 ! This error can invariably be fixed by improving the computation of maxradius 107 ! in the call to init_cube_info 108 ! 109 WRITE (*, *) info%max_rad_ga, radius 110 CPABORT("Called with too large radius.") 111 ENDIF 112 113 ! get the min/max indices of a cube that contains a sphere of the given radius around rp 114 ! if the cell is very non-orthogonal this implies that many useless points are included 115 ! this estimate can be improved (i.e. not box but sphere should be used) 116 lb = HUGE(lb) 117 ub = -HUGE(ub) 118 DO i = -1, 1 119 DO j = -1, 1 120 DO k = -1, 1 121 point(1) = rp(1) + i*radius 122 point(2) = rp(2) + j*radius 123 point(3) = rp(3) + k*radius 124 CALL matvec_3x3(res, info%dh_inv, point) 125 lb = MIN(lb, FLOOR(res)) 126 ub = MAX(ub, CEILING(res)) 127 ENDDO 128 ENDDO 129 ENDDO 130 131 END SUBROUTINE return_cube_nonortho 132 133! ************************************************************************************************** 134!> \brief ... 135!> \param info ... 136!> \param radius ... 137!> \param lb_cube ... 138!> \param ub_cube ... 139!> \param sphere_bounds ... 140! ************************************************************************************************** 141 SUBROUTINE return_cube(info, radius, lb_cube, ub_cube, sphere_bounds) 142 143 TYPE(cube_info_type) :: info 144 REAL(KIND=dp) :: radius 145 INTEGER :: lb_cube(3), ub_cube(3) 146 INTEGER, DIMENSION(:), POINTER :: sphere_bounds 147 148 CHARACTER(LEN=*), PARAMETER :: routineN = 'return_cube', routineP = moduleN//':'//routineN 149 150 INTEGER :: imr 151 152 IF (info%orthorhombic) THEN 153 imr = MAX(1, CEILING(radius/info%drmin)) 154 IF (imr .GT. info%max_radius) THEN 155 ! 156 ! This is an important check. If the required radius for mapping the density is different 157 ! from the actual computed one, (significant) errors can occur. 158 ! This error can invariably be fixed by improving the computation of maxradius 159 ! in the call to init_cube_info 160 ! 161 CPABORT("Called with too large radius.") 162 ENDIF 163 lb_cube(:) = info%lb_cube(:, imr) 164 ub_cube(:) = info%ub_cube(:, imr) 165 sphere_bounds => info%sphere_bounds(imr)%p 166 ELSE 167 ! nothing yet, we should check the radius 168 ENDIF 169 170 END SUBROUTINE return_cube 171 172 ! this is the integer max radius of the cube 173! ************************************************************************************************** 174!> \brief ... 175!> \param info ... 176!> \return ... 177! ************************************************************************************************** 178 INTEGER FUNCTION return_cube_max_iradius(info) 179 TYPE(cube_info_type) :: info 180 181 return_cube_max_iradius = info%max_radius 182 END FUNCTION return_cube_max_iradius 183 184! ************************************************************************************************** 185!> \brief ... 186!> \param info ... 187! ************************************************************************************************** 188 SUBROUTINE destroy_cube_info(info) 189 TYPE(cube_info_type) :: info 190 191 CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_cube_info', & 192 routineP = moduleN//':'//routineN 193 194 INTEGER :: i 195 196 IF (info%orthorhombic) THEN 197 DEALLOCATE (info%lb_cube) 198 DEALLOCATE (info%ub_cube) 199 DEALLOCATE (info%sphere_bounds_count) 200 DO i = 1, info%max_radius 201 DEALLOCATE (info%sphere_bounds(i)%p) 202 END DO 203 DEALLOCATE (info%sphere_bounds) 204 ELSE 205 ! no info to be deallocated 206 ENDIF 207 END SUBROUTINE 208 209! ************************************************************************************************** 210!> \brief ... 211!> \param info ... 212!> \param dr ... 213!> \param dh ... 214!> \param dh_inv ... 215!> \param ortho ... 216!> \param max_radius ... 217! ************************************************************************************************** 218 SUBROUTINE init_cube_info(info, dr, dh, dh_inv, ortho, max_radius) 219 TYPE(cube_info_type), INTENT(OUT) :: info 220 REAL(KIND=dp), INTENT(IN) :: dr(3), dh(3, 3), dh_inv(3, 3) 221 LOGICAL, INTENT(IN) :: ortho 222 REAL(KIND=dp), INTENT(IN) :: max_radius 223 224 CHARACTER(LEN=*), PARAMETER :: routineN = 'init_cube_info', routineP = moduleN//':'//routineN 225 226 INTEGER :: check_1, check_2, handle, i, igmin, imr, & 227 jg, jg2, jgmin, k, kg, kg2, kgmin, & 228 lb(3), ub(3) 229 REAL(KIND=dp) :: drmin, dxi, dy2, dyi, dz2, dzi, radius, & 230 radius2, rp(3) 231 232 CALL timeset(routineN, handle) 233 info%dr = dr 234 info%dh = dh 235 info%dh_inv = dh_inv 236 info%orthorhombic = ortho 237 info%max_rad_ga = max_radius 238 drmin = MINVAL(dr) 239 info%drmin = drmin 240 241 NULLIFY (info%lb_cube, info%ub_cube, & 242 info%sphere_bounds_count, info%sphere_bounds) 243 244 IF (.NOT. info%orthorhombic) THEN 245 246 rp = 0.0_dp 247 ! 248 ! could still be wrong (maybe needs an additional +1 to account for off-gridpoint rp's) 249 ! 250 CALL return_cube_nonortho(info, max_radius, lb, ub, rp) 251 info%max_radius = MAX(MAXVAL(ABS(lb)), MAXVAL(ABS(ub))) 252 253 ELSE 254 255 ! this info is specialized to orthogonal grids 256 imr = CEILING((max_radius)/drmin) 257 info%max_radius = imr 258 dzi = 1.0_dp/dr(3) 259 dyi = 1.0_dp/dr(2) 260 dxi = 1.0_dp/dr(1) 261 dz2 = (dr(3))**2 262 dy2 = (dr(2))**2 263 264 ALLOCATE (info%lb_cube(3, imr), info%ub_cube(3, imr), & 265 info%sphere_bounds_count(imr), info%sphere_bounds(imr)) 266 check_1 = 0 267 check_2 = 0 268! count and allocate 269 270 DO i = 1, imr 271 k = 1 272 radius = i*drmin 273 radius2 = radius**2 274 kgmin = do_and_hide_it_1(dzi, i, drmin, 0.0_dp, 0.0_dp, 0, 0) 275 k = k + 1 276 DO kg = kgmin, 0 277 kg2 = kg*kg 278 jgmin = do_and_hide_it_1(dyi, i, drmin, dz2, 0.0_dp, kg2, 0) 279 k = k + 1 280 DO jg = jgmin, 0 281 jg2 = jg*jg 282 igmin = do_and_hide_it_1(dxi, i, drmin, dz2, dy2, kg2, jg2) 283 check_1 = MODULO((kgmin*97 + jgmin*37 + igmin*113)*check_1 + 1277, 9343) 284 k = k + 1 285 ENDDO 286 ENDDO 287 info%sphere_bounds_count(i) = k - 1 288 ALLOCATE (info%sphere_bounds(i)%p(info%sphere_bounds_count(i))) 289 ENDDO 290 291! init sphere_bounds array 292 ! notice : as many points in lb_cube..0 as 1..ub_cube 293 DO i = 1, imr 294 k = 1 295 radius = i*drmin 296 info%lb_cube(:, i) = -1 297 radius2 = radius**2 298 kgmin = do_and_hide_it_1(dzi, i, drmin, 0.0_dp, 0.0_dp, 0, 0) 299 info%lb_cube(3, i) = MIN(kgmin, info%lb_cube(3, i)) 300 info%sphere_bounds(i)%p(k) = kgmin 301 k = k + 1 302 DO kg = kgmin, 0 303 kg2 = kg*kg 304 jgmin = do_and_hide_it_1(dyi, i, drmin, dz2, 0.0_dp, kg2, 0) 305 info%lb_cube(2, i) = MIN(jgmin, info%lb_cube(2, i)) 306 info%sphere_bounds(i)%p(k) = jgmin 307 k = k + 1 308 DO jg = jgmin, 0 309 jg2 = jg*jg 310 igmin = do_and_hide_it_1(dxi, i, drmin, dz2, dy2, kg2, jg2) 311 check_2 = MODULO((kgmin*97 + jgmin*37 + igmin*113)*check_2 + 1277, 9343) 312 info%lb_cube(1, i) = MIN(igmin, info%lb_cube(1, i)) 313 info%sphere_bounds(i)%p(k) = igmin 314 k = k + 1 315 ENDDO 316 ENDDO 317 info%ub_cube(:, i) = 1 - info%lb_cube(:, i) 318 ENDDO 319 IF (check_1 .NE. check_2) THEN 320 CPABORT("Irreproducible fp math caused memory corruption") 321 END IF 322 323 END IF 324 325 CALL timestop(handle) 326 327 END SUBROUTINE 328 329 ! try to hide things from the optimizer, so that we get the same numbers, 330 ! always (this solves the optimisation problems with the intel and nag compiler 331 ! in which the counting loops and execution loops above are executed a different 332 ! number of times, even at -O1 333! ************************************************************************************************** 334!> \brief ... 335!> \param prefactor ... 336!> \param i ... 337!> \param drmin ... 338!> \param dz2 ... 339!> \param dy2 ... 340!> \param kg2 ... 341!> \param jg2 ... 342!> \return ... 343! ************************************************************************************************** 344 FUNCTION do_and_hide_it_1(prefactor, i, drmin, dz2, dy2, kg2, jg2) RESULT(res) 345 REAL(KIND=dp) :: prefactor 346 INTEGER :: i 347 REAL(KIND=dp) :: drmin, dz2, dy2 348 INTEGER :: kg2, jg2, res 349 350 REAL(KIND=dp), DIMENSION(:), POINTER :: buf 351 352 ALLOCATE (buf(4)) 353 buf(1) = prefactor 354 buf(2) = drmin 355 buf(3) = dz2 356 buf(4) = dy2 357 res = do_and_hide_it_2(buf, i, jg2, kg2) 358 DEALLOCATE (buf) 359 END FUNCTION do_and_hide_it_1 360 361! ************************************************************************************************** 362!> \brief ... 363!> \param buf ... 364!> \param i ... 365!> \param jg2 ... 366!> \param kg2 ... 367!> \return ... 368! ************************************************************************************************** 369 FUNCTION do_and_hide_it_2(buf, i, jg2, kg2) RESULT(res) 370 REAL(KIND=dp), DIMENSION(:), POINTER :: buf 371 INTEGER :: i, jg2, kg2, res 372 373 buf(2) = (i*buf(2))**2 374 res = CEILING(-0.1E-7_dp - buf(1)*SQRT(MAX(buf(2) - kg2*buf(3) - jg2*buf(4), 0.0_dp))) 375 END FUNCTION do_and_hide_it_2 376 377END MODULE cube_utils 378 379