(git:ccc2433)
ps_wavelet_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Definition and initialisation of the ps_wavelet data type.
10 !> \history 01.2014 Renamed from ps_wavelet_types to disentangle dependencies (Ole Schuett)
11 !> \author Florian Schiffmann (09.2007,fschiff)
12 ! **************************************************************************************************
14 
15  USE bibliography, ONLY: genovese2006,&
16  genovese2007,&
17  cite_reference
18  USE kinds, ONLY: dp
20  USE ps_wavelet_types, ONLY: wavelet0d,&
22  ps_wavelet_type
24  psolver,&
27  USE pw_grid_types, ONLY: pw_grid_type
28  USE pw_poisson_types, ONLY: pw_poisson_parameter_type
29  USE pw_types, ONLY: pw_r3d_rs_type
30  USE util, ONLY: get_limit
31 #include "../base/base_uses.f90"
32 
33  IMPLICIT NONE
34 
35  PRIVATE
36 
37  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_wavelet_methods'
38 
39 ! *** Public data types ***
40 
41  PUBLIC :: ps_wavelet_create, &
45 
46 CONTAINS
47 
48 ! **************************************************************************************************
49 !> \brief creates the ps_wavelet_type which is needed for the link to
50 !> the Poisson Solver of Luigi Genovese
51 !> \param poisson_params ...
52 !> \param wavelet wavelet to create
53 !> \param pw_grid the grid that is used to create the wavelet kernel
54 !> \author Flroian Schiffmann
55 ! **************************************************************************************************
56  SUBROUTINE ps_wavelet_create(poisson_params, wavelet, pw_grid)
57  TYPE(pw_poisson_parameter_type), INTENT(IN) :: poisson_params
58  TYPE(ps_wavelet_type), POINTER :: wavelet
59  TYPE(pw_grid_type), POINTER :: pw_grid
60 
61  CHARACTER(len=*), PARAMETER :: routinen = 'ps_wavelet_create'
62 
63  INTEGER :: handle, iproc, nproc, nx, ny, nz
64  REAL(kind=dp) :: hx, hy, hz
65 
66  CALL timeset(routinen, handle)
67 
68  CALL cite_reference(genovese2006)
69  CALL cite_reference(genovese2007)
70 
71  IF (ASSOCIATED(wavelet)) THEN
72  CALL ps_wavelet_release(wavelet)
73  NULLIFY (wavelet)
74  END IF
75 
76  ALLOCATE (wavelet)
77 
78  nx = pw_grid%npts(1)
79  ny = pw_grid%npts(2)
80  nz = pw_grid%npts(3)
81 
82  hx = pw_grid%dr(1)
83  hy = pw_grid%dr(2)
84  hz = pw_grid%dr(3)
85 
86  nproc = product(pw_grid%para%rs_dims)
87 
88  iproc = pw_grid%para%rs_mpo
89 
90  NULLIFY (wavelet%karray, wavelet%rho_z_sliced)
91 
92  wavelet%geocode = poisson_params%wavelet_geocode
93  wavelet%method = poisson_params%wavelet_method
94  wavelet%special_dimension = poisson_params%wavelet_special_dimension
95  wavelet%itype_scf = poisson_params%wavelet_scf_type
96  wavelet%datacode = "D"
97 
98  IF (poisson_params%wavelet_method == wavelet0d) THEN
99  IF (hx .NE. hy) &
100  cpabort("Poisson solver for non cubic cells not yet implemented")
101  IF (hz .NE. hy) &
102  cpabort("Poisson solver for non cubic cells not yet implemented")
103  END IF
104 
105  CALL rs_z_slice_distribution(wavelet, pw_grid)
106 
107  CALL timestop(handle)
108  END SUBROUTINE ps_wavelet_create
109 
110 ! **************************************************************************************************
111 !> \brief ...
112 !> \param wavelet ...
113 !> \param pw_grid ...
114 ! **************************************************************************************************
115  SUBROUTINE rs_z_slice_distribution(wavelet, pw_grid)
116 
117  TYPE(ps_wavelet_type), POINTER :: wavelet
118  TYPE(pw_grid_type), POINTER :: pw_grid
119 
120  CHARACTER(len=*), PARAMETER :: routinen = 'RS_z_slice_distribution'
121 
122  CHARACTER(LEN=1) :: geocode
123  INTEGER :: handle, iproc, m1, m2, m3, md1, md2, &
124  md3, n1, n2, n3, nd1, nd2, nd3, nproc, &
125  nx, ny, nz, z_dim
126  REAL(kind=dp) :: hx, hy, hz
127 
128  CALL timeset(routinen, handle)
129  nproc = product(pw_grid%para%rs_dims)
130  iproc = pw_grid%para%rs_mpo
131  geocode = wavelet%geocode
132  nx = pw_grid%npts(1)
133  ny = pw_grid%npts(2)
134  nz = pw_grid%npts(3)
135  hx = pw_grid%dr(1)
136  hy = pw_grid%dr(2)
137  hz = pw_grid%dr(3)
138 
139  !calculate Dimensions for the z-distributed density and for the kernel
140 
141  IF (geocode == 'P') THEN
142  CALL p_fft_dimensions(nx, ny, nz, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
143  ELSE IF (geocode == 'S') THEN
144  CALL s_fft_dimensions(nx, ny, nz, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
145  ELSE IF (geocode == 'F') THEN
146  CALL f_fft_dimensions(nx, ny, nz, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
147  END IF
148 
149  wavelet%PS_grid(1) = md1
150  wavelet%PS_grid(2) = md3
151  wavelet%PS_grid(3) = md2
152  z_dim = md2/nproc
153  !!!!!!!!! indices y and z are interchanged !!!!!!!
154  ALLOCATE (wavelet%rho_z_sliced(md1, md3, z_dim))
155 
156  CALL createkernel(geocode, nx, ny, nz, hx, hy, hz, wavelet%itype_scf, iproc, nproc, wavelet%karray, &
157  pw_grid%para%rs_group)
158 
159  CALL timestop(handle)
160  END SUBROUTINE rs_z_slice_distribution
161 
162 ! **************************************************************************************************
163 !> \brief ...
164 !> \param density ...
165 !> \param wavelet ...
166 !> \param pw_grid ...
167 ! **************************************************************************************************
168  SUBROUTINE cp2k_distribution_to_z_slices(density, wavelet, pw_grid)
169 
170  TYPE(pw_r3d_rs_type), INTENT(IN) :: density
171  TYPE(ps_wavelet_type), POINTER :: wavelet
172  TYPE(pw_grid_type), POINTER :: pw_grid
173 
174  CHARACTER(len=*), PARAMETER :: routinen = 'cp2k_distribution_to_z_slices'
175 
176  INTEGER :: dest, handle, i, ii, iproc, j, k, l, &
177  local_z_dim, loz, m, m2, md2, nproc, &
178  should_warn
179  INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl, scount, sdispl, tmp
180  INTEGER, DIMENSION(2) :: cart_pos, lox, loy
181  INTEGER, DIMENSION(3) :: lb, ub
182  REAL(kind=dp) :: max_val_low, max_val_up
183  REAL(kind=dp), DIMENSION(:), POINTER :: rbuf, sbuf
184 
185  CALL timeset(routinen, handle)
186 
187  cpassert(ASSOCIATED(wavelet))
188 
189  nproc = product(pw_grid%para%rs_dims)
190  iproc = pw_grid%para%rs_mpo
191  md2 = wavelet%PS_grid(3)
192  m2 = pw_grid%npts(3)
193  lb(:) = pw_grid%bounds_local(1, :)
194  ub(:) = pw_grid%bounds_local(2, :)
195  local_z_dim = max((md2/nproc), 1)
196 
197  ALLOCATE (sbuf(product(pw_grid%npts_local)))
198  ALLOCATE (rbuf(product(wavelet%PS_grid)/nproc))
199  ALLOCATE (scount(nproc), sdispl(nproc), rcount(nproc), rdispl(nproc), tmp(nproc))
200 
201  rbuf = 0.0_dp
202  ii = 1
203  DO k = lb(3), ub(3)
204  DO j = lb(2), ub(2)
205  DO i = lb(1), ub(1)
206  sbuf(ii) = density%array(i, j, k)
207  ii = ii + 1
208  END DO
209  END DO
210  END DO
211 
212  should_warn = 0
213  IF (wavelet%geocode == 'S' .OR. wavelet%geocode == 'F') THEN
214  max_val_low = 0._dp
215  max_val_up = 0._dp
216  IF (lb(2) == pw_grid%bounds(1, 2)) max_val_low = maxval(abs(density%array(:, lb(2), :)))
217  IF (ub(2) == pw_grid%bounds(2, 2)) max_val_up = maxval(abs(density%array(:, ub(2), :)))
218  IF (max_val_low .GE. 0.0001_dp) should_warn = 1
219  IF (max_val_up .GE. 0.0001_dp) should_warn = 1
220  IF (wavelet%geocode == 'F') THEN
221  max_val_low = 0._dp
222  max_val_up = 0._dp
223  IF (lb(1) == pw_grid%bounds(1, 1)) max_val_low = maxval(abs(density%array(lb(1), :, :)))
224  IF (ub(1) == pw_grid%bounds(2, 1)) max_val_up = maxval(abs(density%array(ub(1), :, :)))
225  IF (max_val_low .GE. 0.0001_dp) should_warn = 1
226  IF (max_val_up .GE. 0.0001_dp) should_warn = 1
227  max_val_low = 0._dp
228  max_val_up = 0._dp
229  IF (lb(3) == pw_grid%bounds(1, 3)) max_val_low = maxval(abs(density%array(:, :, lb(3))))
230  IF (ub(3) == pw_grid%bounds(2, 3)) max_val_up = maxval(abs(density%array(:, :, ub(3))))
231  IF (max_val_low .GE. 0.0001_dp) should_warn = 1
232  IF (max_val_up .GE. 0.0001_dp) should_warn = 1
233  END IF
234  END IF
235 
236  CALL pw_grid%para%group%max(should_warn)
237  IF (should_warn > 0 .AND. iproc == 0) &
238  cpwarn("Density non-zero on the edges of the unit cell: wrong results in WAVELET solver")
239 
240  DO i = 0, pw_grid%para%rs_dims(1) - 1
241  DO j = 0, pw_grid%para%rs_dims(2) - 1
242  cart_pos = (/i, j/)
243  CALL pw_grid%para%rs_group%rank_cart(cart_pos, dest)
244  IF ((ub(1) .GE. lb(1)) .AND. (ub(2) .GE. lb(2))) THEN
245  IF (dest*local_z_dim .LE. m2) THEN
246  IF ((dest + 1)*local_z_dim .LE. m2) THEN
247  scount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*local_z_dim)
248  ELSE
249  scount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*mod(m2, local_z_dim))
250  END IF
251  ELSE
252  scount(dest + 1) = 0
253  END IF
254  ELSE
255  scount(dest + 1) = 0
256  END IF
257  lox = get_limit(pw_grid%npts(1), pw_grid%para%rs_dims(1), i)
258  loy = get_limit(pw_grid%npts(2), pw_grid%para%rs_dims(2), j)
259  IF ((lox(2) .GE. lox(1)) .AND. (loy(2) .GE. loy(1))) THEN
260  IF (iproc*local_z_dim .LE. m2) THEN
261  IF ((iproc + 1)*local_z_dim .LE. m2) THEN
262  rcount(dest + 1) = abs((lox(2) - lox(1) + 1)*(loy(2) - loy(1) + 1)*local_z_dim)
263  ELSE
264  rcount(dest + 1) = abs((lox(2) - lox(1) + 1)*(loy(2) - loy(1) + 1)*mod(m2, local_z_dim))
265  END IF
266  ELSE
267  rcount(dest + 1) = 0
268  END IF
269  ELSE
270  rcount(dest + 1) = 0
271  END IF
272 
273  END DO
274  END DO
275  sdispl(1) = 0
276  rdispl(1) = 0
277  DO i = 2, nproc
278  sdispl(i) = sdispl(i - 1) + scount(i - 1)
279  rdispl(i) = rdispl(i - 1) + rcount(i - 1)
280  END DO
281  CALL pw_grid%para%rs_group%alltoall(sbuf, scount, sdispl, rbuf, rcount, rdispl)
282  !!!! and now, how to put the right cubes to the right position!!!!!!
283 
284  wavelet%rho_z_sliced = 0.0_dp
285 
286  DO i = 0, pw_grid%para%rs_dims(1) - 1
287  DO j = 0, pw_grid%para%rs_dims(2) - 1
288  cart_pos = (/i, j/)
289  CALL pw_grid%para%rs_group%rank_cart(cart_pos, dest)
290 
291  lox = get_limit(pw_grid%npts(1), pw_grid%para%rs_dims(1), i)
292  loy = get_limit(pw_grid%npts(2), pw_grid%para%rs_dims(2), j)
293  IF (iproc*local_z_dim .LE. m2) THEN
294  IF ((iproc + 1)*local_z_dim .LE. m2) THEN
295  loz = local_z_dim
296  ELSE
297  loz = mod(m2, local_z_dim)
298  END IF
299  ii = 1
300  DO k = 1, loz
301  DO l = loy(1), loy(2)
302  DO m = lox(1), lox(2)
303  wavelet%rho_z_sliced(m, l, k) = rbuf(ii + rdispl(dest + 1))
304  ii = ii + 1
305  END DO
306  END DO
307  END DO
308  END IF
309  END DO
310  END DO
311 
312  DEALLOCATE (sbuf, rbuf, scount, sdispl, rcount, rdispl, tmp)
313 
314  CALL timestop(handle)
315 
316  END SUBROUTINE cp2k_distribution_to_z_slices
317 
318 ! **************************************************************************************************
319 !> \brief ...
320 !> \param density ...
321 !> \param wavelet ...
322 !> \param pw_grid ...
323 ! **************************************************************************************************
324  SUBROUTINE z_slices_to_cp2k_distribution(density, wavelet, pw_grid)
325 
326  TYPE(pw_r3d_rs_type), INTENT(IN) :: density
327  TYPE(ps_wavelet_type), POINTER :: wavelet
328  TYPE(pw_grid_type), POINTER :: pw_grid
329 
330  INTEGER :: dest, i, ii, iproc, j, k, l, &
331  local_z_dim, loz, m, m2, md2, nproc
332  INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl, scount, sdispl, tmp
333  INTEGER, DIMENSION(2) :: cart_pos, lox, loy, min_x, min_y
334  INTEGER, DIMENSION(3) :: lb, ub
335  REAL(kind=dp), DIMENSION(:), POINTER :: rbuf, sbuf
336 
337  cpassert(ASSOCIATED(wavelet))
338 
339  nproc = product(pw_grid%para%rs_dims)
340  iproc = pw_grid%para%rs_mpo
341  md2 = wavelet%PS_grid(3)
342  m2 = pw_grid%npts(3)
343 
344  lb(:) = pw_grid%bounds_local(1, :)
345  ub(:) = pw_grid%bounds_local(2, :)
346 
347  local_z_dim = max((md2/nproc), 1)
348 
349  ALLOCATE (rbuf(product(pw_grid%npts_local)))
350  ALLOCATE (sbuf(product(wavelet%PS_grid)/nproc))
351  ALLOCATE (scount(nproc), sdispl(nproc), rcount(nproc), rdispl(nproc), tmp(nproc))
352  scount = 0
353  rcount = 0
354  rbuf = 0.0_dp
355  ii = 1
356  IF (iproc*local_z_dim .LE. m2) THEN
357  IF ((iproc + 1)*local_z_dim .LE. m2) THEN
358  loz = local_z_dim
359  ELSE
360  loz = mod(m2, local_z_dim)
361  END IF
362  ELSE
363  loz = 0
364  END IF
365 
366  min_x = get_limit(pw_grid%npts(1), pw_grid%para%rs_dims(1), 0)
367  min_y = get_limit(pw_grid%npts(2), pw_grid%para%rs_dims(2), 0)
368  DO i = 0, pw_grid%para%rs_dims(1) - 1
369  DO j = 0, pw_grid%para%rs_dims(2) - 1
370  cart_pos = (/i, j/)
371  CALL pw_grid%para%rs_group%rank_cart(cart_pos, dest)
372  IF ((ub(1) .GE. lb(1)) .AND. (ub(2) .GE. lb(2))) THEN
373  IF (dest*local_z_dim .LE. m2) THEN
374  IF ((dest + 1)*local_z_dim .LE. m2) THEN
375  rcount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*local_z_dim)
376  ELSE
377  rcount(dest + 1) = abs((ub(1) - lb(1) + 1)*(ub(2) - lb(2) + 1)*mod(m2, local_z_dim))
378  END IF
379  ELSE
380  rcount(dest + 1) = 0
381  END IF
382  ELSE
383  rcount(dest + 1) = 0
384  END IF
385  lox = get_limit(pw_grid%npts(1), pw_grid%para%rs_dims(1), i)
386  loy = get_limit(pw_grid%npts(2), pw_grid%para%rs_dims(2), j)
387  IF ((lox(2) .GE. lox(1)) .AND. (loy(2) .GE. loy(1))) THEN
388  scount(dest + 1) = abs((lox(2) - lox(1) + 1)*(loy(2) - loy(1) + 1)*loz)
389  DO k = lox(1) - min_x(1) + 1, lox(2) - min_x(1) + 1
390  DO l = loy(1) - min_y(1) + 1, loy(2) - min_y(1) + 1
391  DO m = 1, loz
392  sbuf(ii) = wavelet%rho_z_sliced(k, l, m)
393  ii = ii + 1
394  END DO
395  END DO
396  END DO
397  ELSE
398  scount(dest + 1) = 0
399  END IF
400  END DO
401  END DO
402  sdispl(1) = 0
403  rdispl(1) = 0
404  DO i = 2, nproc
405  sdispl(i) = sdispl(i - 1) + scount(i - 1)
406  rdispl(i) = rdispl(i - 1) + rcount(i - 1)
407  END DO
408  CALL pw_grid%para%rs_group%alltoall(sbuf, scount, sdispl, rbuf, rcount, rdispl)
409 
410  !!!! and now, how to put the right cubes to the right position!!!!!!
411 
412  DO i = 0, pw_grid%para%rs_dims(1) - 1
413  DO j = 0, pw_grid%para%rs_dims(2) - 1
414  cart_pos = (/i, j/)
415  CALL pw_grid%para%rs_group%rank_cart(cart_pos, dest)
416  IF (dest*local_z_dim .LE. m2) THEN
417  IF ((dest + 1)*local_z_dim .LE. m2) THEN
418  loz = local_z_dim
419  ELSE
420  loz = mod(m2, local_z_dim)
421  END IF
422  ii = 1
423  IF (lb(3) + (dest*local_z_dim) .LE. ub(3)) THEN
424  DO m = lb(1), ub(1)
425  DO l = lb(2), ub(2)
426  DO k = lb(3) + (dest*local_z_dim), lb(3) + (dest*local_z_dim) + loz - 1
427  density%array(m, l, k) = rbuf(ii + rdispl(dest + 1))
428  ii = ii + 1
429  END DO
430  END DO
431  END DO
432  END IF
433  END IF
434  END DO
435  END DO
436  DEALLOCATE (sbuf, rbuf, scount, sdispl, rcount, rdispl, tmp)
437 
438  END SUBROUTINE z_slices_to_cp2k_distribution
439 
440 ! **************************************************************************************************
441 !> \brief ...
442 !> \param wavelet ...
443 !> \param pw_grid ...
444 ! **************************************************************************************************
445  SUBROUTINE ps_wavelet_solve(wavelet, pw_grid)
446 
447  TYPE(ps_wavelet_type), POINTER :: wavelet
448  TYPE(pw_grid_type), POINTER :: pw_grid
449 
450  CHARACTER(len=*), PARAMETER :: routinen = 'ps_wavelet_solve'
451 
452  CHARACTER(LEN=1) :: geocode
453  INTEGER :: handle, iproc, nproc, nx, ny, nz
454  REAL(kind=dp) :: hx, hy, hz
455 
456  CALL timeset(routinen, handle)
457  nproc = product(pw_grid%para%rs_dims)
458  iproc = pw_grid%para%rs_mpo
459  geocode = wavelet%geocode
460  nx = pw_grid%npts(1)
461  ny = pw_grid%npts(2)
462  nz = pw_grid%npts(3)
463  hx = pw_grid%dr(1)
464  hy = pw_grid%dr(2)
465  hz = pw_grid%dr(3)
466 
467  CALL psolver(geocode, iproc, nproc, nx, ny, nz, hx, hy, hz, &
468  wavelet%rho_z_sliced, wavelet%karray, pw_grid)
469  CALL timestop(handle)
470  END SUBROUTINE ps_wavelet_solve
471 
472 END MODULE ps_wavelet_methods
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public genovese2006
Definition: bibliography.F:43
integer, save, public genovese2007
Definition: bibliography.F:43
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Creates the wavelet kernel for the wavelet based poisson solver.
subroutine, public createkernel(geocode, n01, n02, n03, hx, hy, hz, itype_scf, iproc, nproc, kernel, mpi_group)
Allocate a pointer which corresponds to the zero-padded FFT slice needed for calculating the convolut...
Definition and initialisation of the ps_wavelet data type. \history 01.2014 Renamed from ps_wavelet_t...
subroutine, public ps_wavelet_create(poisson_params, wavelet, pw_grid)
creates the ps_wavelet_type which is needed for the link to the Poisson Solver of Luigi Genovese
subroutine, public ps_wavelet_solve(wavelet, pw_grid)
...
subroutine, public z_slices_to_cp2k_distribution(density, wavelet, pw_grid)
...
subroutine, public cp2k_distribution_to_z_slices(density, wavelet, pw_grid)
...
Definition and initialisation of the ps_wavelet data type.
integer, parameter, public wavelet0d
subroutine, public ps_wavelet_release(wavelet)
...
Performs a wavelet based solution of the Poisson equation.
subroutine, public p_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the convolution for the periodic syste...
subroutine, public psolver(geocode, iproc, nproc, n01, n02, n03, hx, hy, hz, rhopot, karray, pw_grid)
Calculate the Poisson equation $\nabla^2 V(x,y,z)=-4 \pi \rho(x,y,z)$ from a given $\rho$,...
subroutine, public s_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the convolution for the surface system...
subroutine, public f_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the zero-padded convolution.
functions related to the poisson solver on regular grids
All kind of helpful little routines.
Definition: util.F:14
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333