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