Line data Source code
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 numerical operations on real-space grid
10 : !> \par History
11 : !> 12.2014 created [Hossein Bani-Hashemian]
12 : !> \author Mohammad Hossein Bani-Hashemian
13 : ! **************************************************************************************************
14 : MODULE rs_methods
15 :
16 : USE kinds, ONLY: dp
17 : USE pw_grid_types, ONLY: pw_grid_type
18 : USE pw_methods, ONLY: pw_integrate_function,&
19 : pw_scale,&
20 : pw_transfer,&
21 : pw_zero
22 : USE pw_pool_types, ONLY: pw_pool_type
23 : USE pw_types, ONLY: pw_c1d_gs_type,&
24 : pw_r3d_rs_type
25 : USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
26 : realspace_grid_type,&
27 : rs_grid_create,&
28 : rs_grid_release,&
29 : rs_grid_zero,&
30 : transfer_pw2rs,&
31 : transfer_rs2pw
32 : #include "../base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 : PRIVATE
36 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rs_methods'
37 :
38 : PUBLIC derive_fdm_cd3, &
39 : derive_fdm_cd5, &
40 : derive_fdm_cd7, &
41 : setup_grid_axes, &
42 : pw_mollifier
43 :
44 : REAL(dp), PARAMETER, PRIVATE :: small_value = 1.0E-8_dp
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief 2nd order finite difference derivative of a function on realspace grid
50 : !> \param f input function
51 : !> \param df derivative of f
52 : !> \param rs_grid real-space grid
53 : !> \par History:
54 : !> - Creation (15.11.2013,MK)
55 : !> - Refactored and moved here from qs_sccs.F (12.2014, Hossein Bani-Hashemian)
56 : !> \author Matthias Krack (MK)
57 : !> \version 1.0
58 : ! **************************************************************************************************
59 200 : SUBROUTINE derive_fdm_cd3(f, df, rs_grid)
60 :
61 : TYPE(pw_r3d_rs_type), INTENT(IN) :: f
62 : TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: df
63 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid
64 :
65 : CHARACTER(LEN=*), PARAMETER :: routineN = 'derive_fdm_cd3'
66 :
67 : INTEGER :: handle, i, j, k
68 : INTEGER, DIMENSION(3) :: lb, ub
69 : REAL(KIND=dp), DIMENSION(3) :: h
70 400 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: drdx, drdy, drdz, r
71 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
72 7200 : TYPE(realspace_grid_type), DIMENSION(3) :: drs_grid
73 :
74 200 : CALL timeset(routineN, handle)
75 :
76 : ! Setup
77 200 : rs_desc => rs_grid%desc
78 200 : CALL transfer_pw2rs(rs_grid, f)
79 800 : DO i = 1, 3
80 600 : CALL rs_grid_create(drs_grid(i), rs_desc)
81 800 : CALL rs_grid_zero(drs_grid(i))
82 : END DO
83 :
84 800 : lb(1:3) = rs_grid%lb_real(1:3)
85 800 : ub(1:3) = rs_grid%ub_real(1:3)
86 200 : r => rs_grid%r
87 200 : drdx => drs_grid(1)%r
88 200 : drdy => drs_grid(2)%r
89 200 : drdz => drs_grid(3)%r
90 :
91 : ! 3-point stencil central differences
92 800 : h(1:3) = 2.0_dp*f%pw_grid%dr(1:3)
93 : !$OMP PARALLEL DO DEFAULT(NONE) &
94 : !$OMP PRIVATE(i,j,k) &
95 200 : !$OMP SHARED(drdx,drdy,drdz,h,lb,r,ub)
96 : DO k = lb(3), ub(3)
97 : DO j = lb(2), ub(2)
98 : DO i = lb(1), ub(1)
99 : drdx(i, j, k) = (r(i + 1, j, k) - r(i - 1, j, k))/h(1)
100 : drdy(i, j, k) = (r(i, j + 1, k) - r(i, j - 1, k))/h(2)
101 : drdz(i, j, k) = (r(i, j, k + 1) - r(i, j, k - 1))/h(3)
102 : END DO
103 : END DO
104 : END DO
105 : !$OMP END PARALLEL DO
106 :
107 : ! Cleanup
108 800 : DO i = 1, 3
109 600 : CALL transfer_rs2pw(drs_grid(i), df(i))
110 800 : CALL rs_grid_release(drs_grid(i))
111 : END DO
112 :
113 200 : CALL timestop(handle)
114 :
115 400 : END SUBROUTINE derive_fdm_cd3
116 :
117 : ! **************************************************************************************************
118 : !> \brief 4th order finite difference derivative of a function on realspace grid
119 : !> \param f input function
120 : !> \param df derivative of f
121 : !> \param rs_grid real-space grid
122 : !> \par History:
123 : !> - Creation (15.11.2013,MK)
124 : !> - Refactored and moved here from qs_sccs.F (12.2014, Hossein Bani-Hashemian)
125 : !> \author Matthias Krack (MK)
126 : !> \version 1.0
127 : ! **************************************************************************************************
128 192 : SUBROUTINE derive_fdm_cd5(f, df, rs_grid)
129 :
130 : TYPE(pw_r3d_rs_type), INTENT(IN) :: f
131 : TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: df
132 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid
133 :
134 : CHARACTER(LEN=*), PARAMETER :: routineN = 'derive_fdm_cd5'
135 :
136 : INTEGER :: handle, i, j, k
137 : INTEGER, DIMENSION(3) :: lb, ub
138 : REAL(KIND=dp), DIMENSION(3) :: h
139 384 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: drdx, drdy, drdz, r
140 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
141 6912 : TYPE(realspace_grid_type), DIMENSION(3) :: drs_grid
142 :
143 192 : CALL timeset(routineN, handle)
144 :
145 : ! Setup
146 192 : rs_desc => rs_grid%desc
147 192 : CALL transfer_pw2rs(rs_grid, f)
148 768 : DO i = 1, 3
149 576 : CALL rs_grid_create(drs_grid(i), rs_desc)
150 768 : CALL rs_grid_zero(drs_grid(i))
151 : END DO
152 :
153 768 : lb(1:3) = rs_grid%lb_real(1:3)
154 768 : ub(1:3) = rs_grid%ub_real(1:3)
155 192 : r => rs_grid%r
156 192 : drdx => drs_grid(1)%r
157 192 : drdy => drs_grid(2)%r
158 192 : drdz => drs_grid(3)%r
159 :
160 : ! 5-point stencil central differences
161 768 : h(1:3) = 12.0_dp*f%pw_grid%dr(1:3)
162 : !$OMP PARALLEL DO DEFAULT(NONE) &
163 : !$OMP PRIVATE(i,j,k) &
164 192 : !$OMP SHARED(drdx,drdy,drdz,h,lb,r,ub)
165 : DO k = lb(3), ub(3)
166 : DO j = lb(2), ub(2)
167 : DO i = lb(1), ub(1)
168 : drdx(i, j, k) = (r(i - 2, j, k) - r(i + 2, j, k) + 8.0_dp*(r(i + 1, j, k) - r(i - 1, j, k)))/h(1)
169 : drdy(i, j, k) = (r(i, j - 2, k) - r(i, j + 2, k) + 8.0_dp*(r(i, j + 1, k) - r(i, j - 1, k)))/h(2)
170 : drdz(i, j, k) = (r(i, j, k - 2) - r(i, j, k + 2) + 8.0_dp*(r(i, j, k + 1) - r(i, j, k - 1)))/h(3)
171 : END DO
172 : END DO
173 : END DO
174 : !$OMP END PARALLEL DO
175 :
176 : ! Cleanup
177 768 : DO i = 1, 3
178 576 : CALL transfer_rs2pw(drs_grid(i), df(i))
179 768 : CALL rs_grid_release(drs_grid(i))
180 : END DO
181 :
182 192 : CALL timestop(handle)
183 :
184 384 : END SUBROUTINE derive_fdm_cd5
185 :
186 : ! **************************************************************************************************
187 : !> \brief 6th order finite difference derivative of a function on realspace grid
188 : !> \param f input function
189 : !> \param df derivative of f
190 : !> \param rs_grid real-space grid
191 : !> \par History:
192 : !> - Creation (15.11.2013,MK)
193 : !> - Refactored and moved here from qs_sccs.F (12.2014, Hossein Bani-Hashemian)
194 : !> \author Matthias Krack (MK)
195 : !> \version 1.0
196 : ! **************************************************************************************************
197 96 : SUBROUTINE derive_fdm_cd7(f, df, rs_grid)
198 :
199 : TYPE(pw_r3d_rs_type), INTENT(IN) :: f
200 : TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: df
201 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid
202 :
203 : CHARACTER(LEN=*), PARAMETER :: routineN = 'derive_fdm_cd7'
204 :
205 : INTEGER :: handle, i, j, k
206 : INTEGER, DIMENSION(3) :: lb, ub
207 : REAL(KIND=dp), DIMENSION(3) :: h
208 192 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: drdx, drdy, drdz, r
209 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
210 3456 : TYPE(realspace_grid_type), DIMENSION(3) :: drs_grid
211 :
212 96 : CALL timeset(routineN, handle)
213 :
214 : ! Setup
215 96 : rs_desc => rs_grid%desc
216 96 : CALL transfer_pw2rs(rs_grid, f)
217 384 : DO i = 1, 3
218 288 : CALL rs_grid_create(drs_grid(i), rs_desc)
219 384 : CALL rs_grid_zero(drs_grid(i))
220 : END DO
221 :
222 384 : lb(1:3) = rs_grid%lb_real(1:3)
223 384 : ub(1:3) = rs_grid%ub_real(1:3)
224 96 : r => rs_grid%r
225 96 : drdx => drs_grid(1)%r
226 96 : drdy => drs_grid(2)%r
227 96 : drdz => drs_grid(3)%r
228 :
229 : ! 7-point stencil central differences
230 384 : h(1:3) = 60.0_dp*f%pw_grid%dr(1:3)
231 : !$OMP PARALLEL DO DEFAULT(NONE) &
232 : !$OMP PRIVATE(i,j,k) &
233 96 : !$OMP SHARED(drdx,drdy,drdz,h,lb,r,ub)
234 : DO k = lb(3), ub(3)
235 : DO j = lb(2), ub(2)
236 : DO i = lb(1), ub(1)
237 : drdx(i, j, k) = (r(i + 3, j, k) - r(i - 3, j, k) + 9.0_dp*(r(i - 2, j, k) - r(i + 2, j, k)) + &
238 : 45.0_dp*(r(i + 1, j, k) - r(i - 1, j, k)))/h(1)
239 : drdy(i, j, k) = (r(i, j + 3, k) - r(i, j - 3, k) + 9.0_dp*(r(i, j - 2, k) - r(i, j + 2, k)) + &
240 : 45.0_dp*(r(i, j + 1, k) - r(i, j - 1, k)))/h(2)
241 : drdz(i, j, k) = (r(i, j, k + 3) - r(i, j, k - 3) + 9.0_dp*(r(i, j, k - 2) - r(i, j, k + 2)) + &
242 : 45.0_dp*(r(i, j, k + 1) - r(i, j, k - 1)))/h(3)
243 : END DO
244 : END DO
245 : END DO
246 : !$OMP END PARALLEL DO
247 :
248 : ! Cleanup
249 384 : DO i = 1, 3
250 288 : CALL transfer_rs2pw(drs_grid(i), df(i))
251 384 : CALL rs_grid_release(drs_grid(i))
252 : END DO
253 :
254 96 : CALL timestop(handle)
255 :
256 192 : END SUBROUTINE derive_fdm_cd7
257 :
258 : ! **************************************************************************************************
259 : !> \brief returns the global axes and the portion of the axes that are local to
260 : !> the current mpi rank
261 : !> \param pw_grid plane wave grid
262 : !> \param x_glbl x grid vector of the simulation box
263 : !> \param y_glbl y grid vector of the simulation box
264 : !> \param z_glbl z grid vector of the simulation box
265 : !> \param x_locl x grid vector of the simulation box local to this process
266 : !> \param y_locl y grid vector of the simulation box local to this process
267 : !> \param z_locl z grid vector of the simulation box local to this process
268 : !> \par History
269 : !> 07.2014 created [Hossein Bani-Hashemian]
270 : !> 07.2015 moved here from dirichlet_bc_utils.F [Hossein Bani-Hashemian]
271 : !> \author Mohammad Hossein Bani-Hashemian
272 : ! **************************************************************************************************
273 80 : SUBROUTINE setup_grid_axes(pw_grid, x_glbl, y_glbl, z_glbl, x_locl, y_locl, z_locl)
274 :
275 : TYPE(pw_grid_type), INTENT(IN) :: pw_grid
276 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: x_glbl, y_glbl, z_glbl, x_locl, y_locl, &
277 : z_locl
278 :
279 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_grid_axes'
280 :
281 : INTEGER :: glb1, glb2, glb3, gub1, gub2, gub3, &
282 : handle, i, lb1, lb2, lb3, ub1, ub2, ub3
283 80 : INTEGER, ALLOCATABLE, DIMENSION(:) :: gindx, gindy, gindz, lindx, lindy, lindz
284 : INTEGER, DIMENSION(2, 3) :: bounds, bounds_local
285 : INTEGER, DIMENSION(3) :: npts, npts_local
286 : REAL(dp), DIMENSION(3) :: dr
287 :
288 80 : CALL timeset(routineN, handle)
289 :
290 320 : dr = pw_grid%dr
291 800 : bounds = pw_grid%bounds
292 800 : bounds_local = pw_grid%bounds_local
293 320 : npts = pw_grid%npts
294 320 : npts_local = pw_grid%npts_local
295 :
296 : ! local and global lower and upper bounds
297 80 : glb1 = bounds(1, 1); gub1 = bounds(2, 1)
298 80 : glb2 = bounds(1, 2); gub2 = bounds(2, 2)
299 80 : glb3 = bounds(1, 3); gub3 = bounds(2, 3)
300 80 : lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
301 80 : lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
302 80 : lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
303 :
304 560 : ALLOCATE (lindx(lb1:ub1), lindy(lb2:ub2), lindz(lb3:ub3))
305 560 : ALLOCATE (gindx(glb1:gub1), gindy(glb2:gub2), gindz(glb3:gub3))
306 560 : ALLOCATE (x_locl(lb1:ub1), y_locl(lb2:ub2), z_locl(lb3:ub3))
307 560 : ALLOCATE (x_glbl(glb1:gub1), y_glbl(glb2:gub2), z_glbl(glb3:gub3))
308 :
309 10396 : gindx(:) = (/(i, i=0, npts(1) - 1)/)
310 10564 : gindy(:) = (/(i, i=0, npts(2) - 1)/)
311 10340 : gindz(:) = (/(i, i=0, npts(3) - 1)/)
312 5278 : lindx(:) = (/(i, i=0, npts_local(1) - 1)/)
313 10564 : lindy(:) = (/(i, i=0, npts_local(2) - 1)/)
314 10340 : lindz(:) = (/(i, i=0, npts_local(3) - 1)/)
315 :
316 5198 : x_glbl(:) = gindx*dr(1)
317 5282 : y_glbl(:) = gindy*dr(2)
318 5170 : z_glbl(:) = gindz*dr(3)
319 :
320 : ! map [0 .. (npts_local-1)] --> [lb .. ub]
321 80 : IF (lb1 .EQ. ub1) THEN
322 0 : lindx(:) = lb1
323 : ELSE
324 2639 : lindx(:) = lindx(:)*((ub1 - lb1)/(npts_local(1) - 1)) + lb1
325 : END IF
326 80 : IF (lb2 .EQ. ub2) THEN
327 0 : lindy(:) = lb2
328 : ELSE
329 5282 : lindy(:) = lindy(:)*((ub2 - lb2)/(npts_local(2) - 1)) + lb2
330 : END IF
331 80 : IF (lb3 .EQ. ub3) THEN
332 0 : lindz(:) = lb3
333 : ELSE
334 5170 : lindz(:) = lindz(:)*((ub3 - lb3)/(npts_local(3) - 1)) + lb3
335 : END IF
336 :
337 2639 : x_locl(:) = x_glbl(lindx)
338 5282 : y_locl(:) = y_glbl(lindy)
339 5170 : z_locl(:) = z_glbl(lindz)
340 :
341 80 : CALL timestop(handle)
342 :
343 160 : END SUBROUTINE setup_grid_axes
344 :
345 : ! **************************************************************************************************
346 : !> \brief convolutes a function with a smoothing kernel K_\zeta
347 : !> v * K_\zeta
348 : !> K_\zeta is the standard mollifier defined as:
349 : !> K_\zeta(x) = \frac{1}{\zeta^3} K(\frac{x}{\zeta})
350 : !> where
351 : !> K(x) = \kappa \exp (\frac{1}{|x|^2 - 1}), if |x| <= 1
352 : !> = 0, otherwise
353 : !> \param pw_pool pool of pw grid
354 : !> \param zeta parameter \zeta defining the width of the mollifier
355 : !> \param x_glbl x grid vector of the simulation box
356 : !> \param y_glbl y grid vector of the simulation box
357 : !> \param z_glbl z grid vector of the simulation box
358 : !> \param pw_in the input function
359 : !> \param pw_out the convoluted function
360 : !> \par History
361 : !> 10.2014 created [Hossein Bani-Hashemian]
362 : !> 07.2015 moved here from ps_implicit_methods.F [Hossein Bani-Hashemian]
363 : !> \author Mohammad Hossein Bani-Hashemian
364 : ! **************************************************************************************************
365 28 : SUBROUTINE pw_mollifier(pw_pool, zeta, x_glbl, y_glbl, z_glbl, pw_in, pw_out)
366 :
367 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
368 : REAL(dp), INTENT(IN) :: zeta
369 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_glbl, y_glbl, z_glbl
370 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
371 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
372 :
373 : CHARACTER(LEN=*), PARAMETER :: routineN = 'pw_mollifier'
374 :
375 : INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
376 : ub2, ub3
377 : INTEGER, DIMENSION(2, 3) :: bounds, bounds_local
378 : REAL(dp) :: normfact, xi, xmax, xmin, yj, ymax, &
379 : ymin, zk, zmax, zmin
380 : REAL(dp), DIMENSION(3, 3) :: dh
381 : TYPE(pw_c1d_gs_type) :: G_gs, pw_in_gs, pw_out_gs
382 : TYPE(pw_grid_type), POINTER :: pw_grid
383 : TYPE(pw_r3d_rs_type) :: G
384 :
385 28 : CALL timeset(routineN, handle)
386 :
387 28 : pw_grid => pw_in%pw_grid
388 364 : dh = pw_grid%dh
389 280 : bounds_local = pw_grid%bounds_local
390 280 : bounds = pw_grid%bounds
391 :
392 28 : lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
393 28 : lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
394 28 : lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
395 :
396 28 : CALL pw_pool%create_pw(G)
397 28 : CALL pw_pool%create_pw(G_gs)
398 28 : CALL pw_pool%create_pw(pw_in_gs)
399 28 : CALL pw_pool%create_pw(pw_out_gs)
400 :
401 28 : CALL pw_zero(G)
402 28 : xmin = x_glbl(bounds(1, 1)); xmax = x_glbl(bounds(2, 1))
403 28 : ymin = y_glbl(bounds(1, 2)); ymax = y_glbl(bounds(2, 2))
404 28 : zmin = z_glbl(bounds(1, 3)); zmax = z_glbl(bounds(2, 3))
405 :
406 2044 : DO k = lb3, ub3
407 147196 : DO j = lb2, ub2
408 5372640 : DO i = lb1, ub1
409 5225472 : xi = x_glbl(i); yj = y_glbl(j); zk = z_glbl(k)
410 21047040 : IF (norm2((/(xi - xmin), (yj - ymin), (zk - zmin)/)) .LT. zeta - small_value) THEN
411 5940 : G%array(i, j, k) = EXP(1.0_dp/(norm2((/(xi - xmin), (yj - ymin), (zk - zmin)/)/zeta)**2 - 1))
412 20895948 : ELSE IF (norm2((/(xi - xmax), (yj - ymax), (zk - zmax)/)) .LT. zeta - small_value) THEN
413 5940 : G%array(i, j, k) = EXP(1.0_dp/(norm2((/(xi - xmax), (yj - ymax), (zk - zmax)/)/zeta)**2 - 1))
414 20890008 : ELSE IF (norm2((/(xi - xmin), (yj - ymax), (zk - zmax)/)) .LT. zeta - small_value) THEN
415 5940 : G%array(i, j, k) = EXP(1.0_dp/(norm2((/(xi - xmin), (yj - ymax), (zk - zmax)/)/zeta)**2 - 1))
416 20884068 : ELSE IF (norm2((/(xi - xmax), (yj - ymin), (zk - zmax)/)) .LT. zeta - small_value) THEN
417 5940 : G%array(i, j, k) = EXP(1.0_dp/(norm2((/(xi - xmax), (yj - ymin), (zk - zmax)/)/zeta)**2 - 1))
418 20878128 : ELSE IF (norm2((/(xi - xmax), (yj - ymax), (zk - zmin)/)) .LT. zeta - small_value) THEN
419 5940 : G%array(i, j, k) = EXP(1.0_dp/(norm2((/(xi - xmax), (yj - ymax), (zk - zmin)/)/zeta)**2 - 1))
420 20872188 : ELSE IF (norm2((/(xi - xmin), (yj - ymin), (zk - zmax)/)) .LT. zeta - small_value) THEN
421 5940 : G%array(i, j, k) = EXP(1.0_dp/(norm2((/(xi - xmin), (yj - ymin), (zk - zmax)/)/zeta)**2 - 1))
422 20866248 : ELSE IF (norm2((/(xi - xmin), (yj - ymax), (zk - zmin)/)) .LT. zeta - small_value) THEN
423 5940 : G%array(i, j, k) = EXP(1.0_dp/(norm2((/(xi - xmin), (yj - ymax), (zk - zmin)/)/zeta)**2 - 1))
424 20860308 : ELSE IF (norm2((/(xi - xmax), (yj - ymin), (zk - zmin)/)) .LT. zeta - small_value) THEN
425 5940 : G%array(i, j, k) = EXP(1.0_dp/(norm2((/(xi - xmax), (yj - ymin), (zk - zmin)/)/zeta)**2 - 1))
426 : END IF
427 : END DO
428 : END DO
429 : END DO
430 28 : CALL pw_scale(G, (1.0_dp/zeta)**3)
431 28 : normfact = pw_integrate_function(G)
432 28 : CALL pw_scale(G, 1.0_dp/normfact)
433 :
434 28 : CALL pw_transfer(G, G_gs)
435 28 : CALL pw_transfer(pw_in, pw_in_gs)
436 5225500 : pw_out_gs%array = G_gs%array*pw_in_gs%array
437 28 : CALL pw_transfer(pw_out_gs, pw_out)
438 :
439 : ! multiply by the reciprocal of the forward Fourier transform normalization prefactor (here 1/N, by convention)
440 28 : CALL pw_scale(pw_out, REAL(pw_grid%ngpts, KIND=dp))
441 : ! from discrete convolution to continuous convolution
442 28 : CALL pw_scale(pw_out, pw_grid%dvol)
443 :
444 2044 : DO k = lb3, ub3
445 147196 : DO j = lb2, ub2
446 5372640 : DO i = lb1, ub1
447 5370624 : IF (ABS(pw_out%array(i, j, k)) .LE. 1.0E-10_dp) pw_out%array(i, j, k) = 0.0_dp
448 : END DO
449 : END DO
450 : END DO
451 :
452 28 : CALL pw_pool%give_back_pw(G)
453 28 : CALL pw_pool%give_back_pw(G_gs)
454 28 : CALL pw_pool%give_back_pw(pw_in_gs)
455 28 : CALL pw_pool%give_back_pw(pw_out_gs)
456 28 : CALL timestop(handle)
457 :
458 28 : END SUBROUTINE pw_mollifier
459 :
460 : END MODULE rs_methods
|
