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 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 7908156 : 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 7908156 : zetp = zeta + zetb
77 31632624 : rp(:) = ra(:) + zetb/zetp*rab(:)
78 134438652 : cube_center(:) = FLOOR(MATMUL(rs_desc%dh_inv, rp))
79 :
80 7908156 : 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 440561 : 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 440561 : 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 0 : WRITE (*, *) info%max_rad_ga, radius
108 0 : CPABORT("Called with too large radius.")
109 : END IF
110 :
111 : ! get the min/max indices of a cube that contains a sphere of the given radius around rp
112 : ! if the cell is very non-orthogonal this implies that many useless points are included
113 : ! this estimate can be improved (i.e. not box but sphere should be used)
114 1762244 : lb = HUGE(lb)
115 1762244 : ub = -HUGE(ub)
116 1762244 : DO i = -1, 1
117 5727293 : DO j = -1, 1
118 17181879 : DO k = -1, 1
119 11895147 : point(1) = rp(1) + i*radius
120 11895147 : point(2) = rp(2) + j*radius
121 11895147 : point(3) = rp(3) + k*radius
122 154636911 : res = MATMUL(info%dh_inv, point)
123 47580588 : lb = MIN(lb, FLOOR(res))
124 51545637 : ub = MAX(ub, CEILING(res))
125 : END DO
126 : END DO
127 : END DO
128 :
129 440561 : END SUBROUTINE return_cube_nonortho
130 :
131 : ! **************************************************************************************************
132 : !> \brief ...
133 : !> \param info ...
134 : !> \param radius ...
135 : !> \param lb_cube ...
136 : !> \param ub_cube ...
137 : !> \param sphere_bounds ...
138 : ! **************************************************************************************************
139 7507840 : SUBROUTINE return_cube(info, radius, lb_cube, ub_cube, sphere_bounds)
140 :
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
145 :
146 : INTEGER :: imr
147 :
148 7507840 : IF (info%orthorhombic) THEN
149 7507840 : imr = MAX(1, CEILING(radius/info%drmin))
150 7507840 : IF (imr .GT. info%max_radius) THEN
151 : !
152 : ! This is an important check. If the required radius for mapping the density is different
153 : ! from the actual computed one, (significant) errors can occur.
154 : ! This error can invariably be fixed by improving the computation of maxradius
155 : ! in the call to init_cube_info
156 : !
157 0 : CPABORT("Called with too large radius.")
158 : END IF
159 30031360 : lb_cube(:) = info%lb_cube(:, imr)
160 30031360 : ub_cube(:) = info%ub_cube(:, imr)
161 7507840 : sphere_bounds => info%sphere_bounds(imr)%p
162 : ELSE
163 : ! nothing yet, we should check the radius
164 : END IF
165 :
166 7507840 : END SUBROUTINE return_cube
167 :
168 : ! this is the integer max radius of the cube
169 : ! **************************************************************************************************
170 : !> \brief ...
171 : !> \param info ...
172 : !> \return ...
173 : ! **************************************************************************************************
174 29366 : INTEGER FUNCTION return_cube_max_iradius(info)
175 : TYPE(cube_info_type) :: info
176 :
177 29366 : return_cube_max_iradius = info%max_radius
178 29366 : END FUNCTION return_cube_max_iradius
179 :
180 : ! **************************************************************************************************
181 : !> \brief ...
182 : !> \param info ...
183 : ! **************************************************************************************************
184 30326 : SUBROUTINE destroy_cube_info(info)
185 : TYPE(cube_info_type) :: info
186 :
187 : INTEGER :: i
188 :
189 30326 : IF (info%orthorhombic) THEN
190 28146 : DEALLOCATE (info%lb_cube)
191 28146 : DEALLOCATE (info%ub_cube)
192 28146 : DEALLOCATE (info%sphere_bounds_count)
193 462064 : DO i = 1, info%max_radius
194 462064 : DEALLOCATE (info%sphere_bounds(i)%p)
195 : END DO
196 28146 : DEALLOCATE (info%sphere_bounds)
197 : ELSE
198 : ! no info to be deallocated
199 : END IF
200 30326 : END SUBROUTINE
201 :
202 : ! **************************************************************************************************
203 : !> \brief ...
204 : !> \param info ...
205 : !> \param dr ...
206 : !> \param dh ...
207 : !> \param dh_inv ...
208 : !> \param ortho ...
209 : !> \param max_radius ...
210 : ! **************************************************************************************************
211 909780 : SUBROUTINE init_cube_info(info, dr, dh, dh_inv, ortho, max_radius)
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
216 :
217 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_cube_info'
218 :
219 : INTEGER :: check_1, check_2, handle, i, igmin, imr, &
220 : jg, jg2, jgmin, k, kg, kg2, kgmin, &
221 : lb(3), ub(3)
222 : REAL(KIND=dp) :: drmin, dxi, dy2, dyi, dz2, dzi, radius, &
223 : radius2, rp(3)
224 :
225 30326 : CALL timeset(routineN, handle)
226 121304 : info%dr = dr
227 394238 : info%dh = dh
228 394238 : info%dh_inv = dh_inv
229 30326 : info%orthorhombic = ortho
230 30326 : info%max_rad_ga = max_radius
231 151630 : drmin = MINVAL(dr)
232 30326 : info%drmin = drmin
233 :
234 30326 : NULLIFY (info%lb_cube, info%ub_cube, &
235 30326 : info%sphere_bounds_count, info%sphere_bounds)
236 :
237 30326 : IF (.NOT. info%orthorhombic) THEN
238 :
239 2180 : rp = 0.0_dp
240 : !
241 : ! could still be wrong (maybe needs an additional +1 to account for off-gridpoint rp's)
242 : !
243 2180 : CALL return_cube_nonortho(info, max_radius, lb, ub, rp)
244 17440 : info%max_radius = MAX(MAXVAL(ABS(lb)), MAXVAL(ABS(ub)))
245 :
246 : ELSE
247 :
248 : ! this info is specialized to orthogonal grids
249 28146 : imr = CEILING((max_radius)/drmin)
250 28146 : info%max_radius = imr
251 28146 : dzi = 1.0_dp/dr(3)
252 28146 : dyi = 1.0_dp/dr(2)
253 28146 : dxi = 1.0_dp/dr(1)
254 28146 : dz2 = (dr(3))**2
255 28146 : dy2 = (dr(2))**2
256 :
257 : ALLOCATE (info%lb_cube(3, imr), info%ub_cube(3, imr), &
258 685952 : info%sphere_bounds_count(imr), info%sphere_bounds(imr))
259 28146 : check_1 = 0
260 28146 : check_2 = 0
261 : ! count and allocate
262 :
263 462064 : DO i = 1, imr
264 433918 : k = 1
265 433918 : radius = i*drmin
266 433918 : radius2 = radius**2
267 433918 : kgmin = do_and_hide_it_1(dzi, i, drmin, 0.0_dp, 0.0_dp, 0, 0)
268 433918 : k = k + 1
269 6647756 : DO kg = kgmin, 0
270 6213838 : kg2 = kg*kg
271 6213838 : jgmin = do_and_hide_it_1(dyi, i, drmin, dz2, 0.0_dp, kg2, 0)
272 6213838 : k = k + 1
273 286838474 : DO jg = jgmin, 0
274 280190718 : jg2 = jg*jg
275 280190718 : igmin = do_and_hide_it_1(dxi, i, drmin, dz2, dy2, kg2, jg2)
276 280190718 : check_1 = MODULO((kgmin*97 + jgmin*37 + igmin*113)*check_1 + 1277, 9343)
277 286404556 : k = k + 1
278 : END DO
279 : END DO
280 433918 : info%sphere_bounds_count(i) = k - 1
281 1329900 : ALLOCATE (info%sphere_bounds(i)%p(info%sphere_bounds_count(i)))
282 : END DO
283 :
284 : ! init sphere_bounds array
285 : ! notice : as many points in lb_cube..0 as 1..ub_cube
286 462064 : DO i = 1, imr
287 433918 : k = 1
288 433918 : radius = i*drmin
289 1735672 : info%lb_cube(:, i) = -1
290 433918 : radius2 = radius**2
291 433918 : kgmin = do_and_hide_it_1(dzi, i, drmin, 0.0_dp, 0.0_dp, 0, 0)
292 433918 : info%lb_cube(3, i) = MIN(kgmin, info%lb_cube(3, i))
293 433918 : info%sphere_bounds(i)%p(k) = kgmin
294 433918 : k = k + 1
295 6647756 : DO kg = kgmin, 0
296 6213838 : kg2 = kg*kg
297 6213838 : jgmin = do_and_hide_it_1(dyi, i, drmin, dz2, 0.0_dp, kg2, 0)
298 6213838 : info%lb_cube(2, i) = MIN(jgmin, info%lb_cube(2, i))
299 6213838 : info%sphere_bounds(i)%p(k) = jgmin
300 6213838 : k = k + 1
301 286838474 : DO jg = jgmin, 0
302 280190718 : jg2 = jg*jg
303 280190718 : igmin = do_and_hide_it_1(dxi, i, drmin, dz2, dy2, kg2, jg2)
304 280190718 : check_2 = MODULO((kgmin*97 + jgmin*37 + igmin*113)*check_2 + 1277, 9343)
305 280190718 : info%lb_cube(1, i) = MIN(igmin, info%lb_cube(1, i))
306 280190718 : info%sphere_bounds(i)%p(k) = igmin
307 286404556 : k = k + 1
308 : END DO
309 : END DO
310 1763818 : info%ub_cube(:, i) = 1 - info%lb_cube(:, i)
311 : END DO
312 28146 : IF (check_1 .NE. check_2) THEN
313 0 : CPABORT("Irreproducible fp math caused memory corruption")
314 : END IF
315 :
316 : END IF
317 :
318 30326 : CALL timestop(handle)
319 :
320 30326 : END SUBROUTINE
321 :
322 : ! try to hide things from the optimizer, so that we get the same numbers,
323 : ! always (this solves the optimisation problems with the intel and nag compiler
324 : ! in which the counting loops and execution loops above are executed a different
325 : ! number of times, even at -O1
326 : ! **************************************************************************************************
327 : !> \brief ...
328 : !> \param prefactor ...
329 : !> \param i ...
330 : !> \param drmin ...
331 : !> \param dz2 ...
332 : !> \param dy2 ...
333 : !> \param kg2 ...
334 : !> \param jg2 ...
335 : !> \return ...
336 : ! **************************************************************************************************
337 573676948 : FUNCTION do_and_hide_it_1(prefactor, i, drmin, dz2, dy2, kg2, jg2) RESULT(res)
338 : REAL(KIND=dp) :: prefactor
339 : INTEGER :: i
340 : REAL(KIND=dp) :: drmin, dz2, dy2
341 : INTEGER :: kg2, jg2, res
342 :
343 573676948 : REAL(KIND=dp), DIMENSION(:), POINTER :: buf
344 :
345 573676948 : ALLOCATE (buf(4))
346 573676948 : buf(1) = prefactor
347 573676948 : buf(2) = drmin
348 573676948 : buf(3) = dz2
349 573676948 : buf(4) = dy2
350 573676948 : res = do_and_hide_it_2(buf, i, jg2, kg2)
351 573676948 : DEALLOCATE (buf)
352 573676948 : END FUNCTION do_and_hide_it_1
353 :
354 : ! **************************************************************************************************
355 : !> \brief ...
356 : !> \param buf ...
357 : !> \param i ...
358 : !> \param jg2 ...
359 : !> \param kg2 ...
360 : !> \return ...
361 : ! **************************************************************************************************
362 573676948 : 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
365 :
366 573676948 : buf(2) = (i*buf(2))**2
367 573676948 : res = CEILING(-0.1E-7_dp - buf(1)*SQRT(MAX(buf(2) - kg2*buf(3) - jg2*buf(4), 0.0_dp)))
368 573676948 : END FUNCTION do_and_hide_it_2
369 :
370 0 : END MODULE cube_utils
371 :
|
