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