(git:374b731)
Loading...
Searching...
No Matches
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,&
17 cite_reference
18 USE kinds, ONLY: dp
20 USE ps_wavelet_types, ONLY: wavelet0d,&
24 psolver,&
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
46CONTAINS
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
472END MODULE ps_wavelet_methods
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public genovese2006
integer, save, public genovese2007
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
parameters for the poisson solver independet of input_section