22 #include "../base/base_uses.f90"
26 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'cube_utils'
28 PUBLIC :: cube_info_type
35 INTEGER,
POINTER,
DIMENSION(:) :: p => null()
38 TYPE :: cube_info_type
39 INTEGER :: max_radius = 0.0_dp
40 REAL(KIND=
dp) :: dr(3) = 0.0_dp, drmin = 0.0_dp
41 REAL(KIND=
dp) :: dh(3, 3) = 0.0_dp
42 REAL(KIND=
dp) :: dh_inv(3, 3) = 0.0_dp
43 LOGICAL :: orthorhombic = .true.
44 INTEGER,
POINTER :: lb_cube(:, :) => null()
45 INTEGER,
POINTER :: ub_cube(:, :) => null()
46 TYPE(cube_ptr),
POINTER,
DIMENSION(:) :: sphere_bounds => null()
47 INTEGER,
POINTER :: sphere_bounds_count(:) => null()
48 REAL(KIND=
dp) :: max_rad_ga = 0.0_dp
49 END TYPE cube_info_type
69 INTEGER,
DIMENSION(3),
INTENT(OUT) :: cube_center
70 TYPE(realspace_grid_desc_type),
POINTER :: rs_desc
71 REAL(kind=
dp),
INTENT(IN) :: zeta, zetb, ra(3), rab(3)
74 REAL(kind=
dp),
DIMENSION(3) :: rp
77 rp(:) = ra(:) + zetb/zetp*rab(:)
78 cube_center(:) = floor(matmul(rs_desc%dh_inv, rp))
92 TYPE(cube_info_type),
INTENT(IN) :: info
93 REAL(kind=
dp),
INTENT(IN) :: radius
94 INTEGER,
INTENT(OUT) :: lb(3), ub(3)
95 REAL(kind=
dp),
INTENT(IN) :: rp(3)
98 REAL(kind=
dp) :: point(3), res(3)
100 IF (radius > info%max_rad_ga)
THEN
107 WRITE (*, *) info%max_rad_ga, radius
108 cpabort(
"Called with too large radius.")
119 point(1) = rp(1) + i*radius
120 point(2) = rp(2) + j*radius
121 point(3) = rp(3) + k*radius
122 res = matmul(info%dh_inv, point)
123 lb = min(lb, floor(res))
124 ub = max(ub, ceiling(res))
139 SUBROUTINE return_cube(info, radius, lb_cube, ub_cube, sphere_bounds)
141 TYPE(cube_info_type) :: info
142 REAL(kind=
dp) :: radius
143 INTEGER :: lb_cube(3), ub_cube(3)
144 INTEGER,
DIMENSION(:),
POINTER :: sphere_bounds
148 IF (info%orthorhombic)
THEN
149 imr = max(1, ceiling(radius/info%drmin))
150 IF (imr .GT. info%max_radius)
THEN
157 cpabort(
"Called with too large radius.")
159 lb_cube(:) = info%lb_cube(:, imr)
160 ub_cube(:) = info%ub_cube(:, imr)
161 sphere_bounds => info%sphere_bounds(imr)%p
175 TYPE(cube_info_type) :: info
185 TYPE(cube_info_type) :: info
189 IF (info%orthorhombic)
THEN
190 DEALLOCATE (info%lb_cube)
191 DEALLOCATE (info%ub_cube)
192 DEALLOCATE (info%sphere_bounds_count)
193 DO i = 1, info%max_radius
194 DEALLOCATE (info%sphere_bounds(i)%p)
196 DEALLOCATE (info%sphere_bounds)
212 TYPE(cube_info_type),
INTENT(OUT) :: info
213 REAL(kind=
dp),
INTENT(IN) :: dr(3), dh(3, 3), dh_inv(3, 3)
214 LOGICAL,
INTENT(IN) :: ortho
215 REAL(kind=
dp),
INTENT(IN) :: max_radius
217 CHARACTER(LEN=*),
PARAMETER :: routinen =
'init_cube_info'
219 INTEGER :: check_1, check_2, handle, i, igmin, imr, &
220 jg, jg2, jgmin, k, kg, kg2, kgmin, &
222 REAL(kind=
dp) :: drmin, dxi, dy2, dyi, dz2, dzi, radius, &
225 CALL timeset(routinen, handle)
229 info%orthorhombic = ortho
230 info%max_rad_ga = max_radius
234 NULLIFY (info%lb_cube, info%ub_cube, &
235 info%sphere_bounds_count, info%sphere_bounds)
237 IF (.NOT. info%orthorhombic)
THEN
244 info%max_radius = max(maxval(abs(lb)), maxval(abs(ub)))
249 imr = ceiling((max_radius)/drmin)
250 info%max_radius = imr
257 ALLOCATE (info%lb_cube(3, imr), info%ub_cube(3, imr), &
258 info%sphere_bounds_count(imr), info%sphere_bounds(imr))
267 kgmin = do_and_hide_it_1(dzi, i, drmin, 0.0_dp, 0.0_dp, 0, 0)
271 jgmin = do_and_hide_it_1(dyi, i, drmin, dz2, 0.0_dp, kg2, 0)
275 igmin = do_and_hide_it_1(dxi, i, drmin, dz2, dy2, kg2, jg2)
276 check_1 =
modulo((kgmin*97 + jgmin*37 + igmin*113)*check_1 + 1277, 9343)
280 info%sphere_bounds_count(i) = k - 1
281 ALLOCATE (info%sphere_bounds(i)%p(info%sphere_bounds_count(i)))
289 info%lb_cube(:, i) = -1
291 kgmin = do_and_hide_it_1(dzi, i, drmin, 0.0_dp, 0.0_dp, 0, 0)
292 info%lb_cube(3, i) = min(kgmin, info%lb_cube(3, i))
293 info%sphere_bounds(i)%p(k) = kgmin
297 jgmin = do_and_hide_it_1(dyi, i, drmin, dz2, 0.0_dp, kg2, 0)
298 info%lb_cube(2, i) = min(jgmin, info%lb_cube(2, i))
299 info%sphere_bounds(i)%p(k) = jgmin
303 igmin = do_and_hide_it_1(dxi, i, drmin, dz2, dy2, kg2, jg2)
304 check_2 =
modulo((kgmin*97 + jgmin*37 + igmin*113)*check_2 + 1277, 9343)
305 info%lb_cube(1, i) = min(igmin, info%lb_cube(1, i))
306 info%sphere_bounds(i)%p(k) = igmin
310 info%ub_cube(:, i) = 1 - info%lb_cube(:, i)
312 IF (check_1 .NE. check_2)
THEN
313 cpabort(
"Irreproducible fp math caused memory corruption")
318 CALL timestop(handle)
337 FUNCTION do_and_hide_it_1(prefactor, i, drmin, dz2, dy2, kg2, jg2)
RESULT(res)
338 REAL(kind=
dp) :: prefactor
340 REAL(kind=
dp) :: drmin, dz2, dy2
341 INTEGER :: kg2, jg2, res
343 REAL(kind=
dp),
DIMENSION(:),
POINTER :: buf
350 res = do_and_hide_it_2(buf, i, jg2, kg2)
352 END FUNCTION do_and_hide_it_1
362 FUNCTION do_and_hide_it_2(buf, i, jg2, kg2)
RESULT(res)
363 REAL(kind=
dp),
DIMENSION(:),
POINTER :: buf
364 INTEGER :: i, jg2, kg2, res
366 buf(2) = (i*buf(2))**2
367 res = ceiling(-0.1e-7_dp - buf(1)*sqrt(max(buf(2) - kg2*buf(3) - jg2*buf(4), 0.0_dp)))
368 END FUNCTION do_and_hide_it_2
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
for a given dr()/dh(r) this will provide the bounds to be used if one wants to go over a sphere-subre...
subroutine, public compute_cube_center(cube_center, rs_desc, zeta, zetb, ra, rab)
unifies the computation of the cube center, so that differences in implementation,...
subroutine, public return_cube(info, radius, lb_cube, ub_cube, sphere_bounds)
...
integer function, public return_cube_max_iradius(info)
...
subroutine, public destroy_cube_info(info)
...
subroutine, public init_cube_info(info, dr, dh, dh_inv, ortho, max_radius)
...
subroutine, public return_cube_nonortho(info, radius, lb, ub, rp)
...
Defines the basic variable types.
integer, parameter, public dp