Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 2 : PROGRAM pw_unittest
9 2 : USE cp_log_handling, ONLY: cp_get_default_logger,&
10 : cp_logger_type
11 : USE kinds, ONLY: dp
12 : USE libcp2k, ONLY: cp2k_finalize,&
13 : cp2k_init
14 : USE mathconstants, ONLY: fourpi,&
15 : pi
16 : USE message_passing, ONLY: mp_comm_type,&
17 : mp_comm_world,&
18 : mp_para_env_type
19 : USE mt_util, ONLY: mt2d_kernel_strain
20 : USE pw_copy_all, ONLY: pw_copy_match
21 : USE pw_grid_types, ONLY: FULLSPACE,&
22 : HALFSPACE,&
23 : pw_grid_type
24 : USE pw_grids, ONLY: pw_grid_create,&
25 : pw_grid_release
26 : USE pw_methods, ONLY: do_accurate_sum,&
27 : do_standard_sum,&
28 : pw_integral_ab
29 : USE pw_poisson_methods, ONLY: pw_poisson_set,&
30 : pw_poisson_solve
31 : USE pw_poisson_types, ONLY: pw_poisson_mt,&
32 : pw_poisson_type
33 : USE pw_pool_types, ONLY: pw_pool_create
34 : USE pw_types, ONLY: pw_c1d_gs_type,&
35 : pw_r1d_gs_type
36 : USE realspace_grid_cube, ONLY: cube_read_values,&
37 : cube_value_format,&
38 : cube_values_format
39 :
40 : IMPLICIT NONE
41 :
42 : INTEGER :: mode
43 : REAL(dp) :: max_error
44 : TYPE(cp_logger_type), POINTER :: logger
45 : TYPE(mp_para_env_type), POINTER :: para_env
46 :
47 2 : CALL cp2k_init()
48 2 : logger => cp_get_default_logger()
49 2 : para_env => logger%para_env
50 2 : CALL check_cube_values()
51 20 : DO mode = 0, 8
52 20 : CALL check_copy(mode)
53 : END DO
54 2 : max_error = 0.0_dp
55 2 : CALL check_pw_integral_and_poisson()
56 2 : CALL check_mt2d_kernel()
57 2 : PRINT '(A,ES16.7)', 'PW integral and MT2D kernel maximum scaled error: ', max_error
58 2 : CALL cp2k_finalize()
59 :
60 : CONTAINS
61 :
62 : ! **************************************************************************************************
63 : !> \brief Preserve three-digit cube exponents and parse adjacent values across z-slices.
64 : ! **************************************************************************************************
65 2 : SUBROUTINE check_cube_values()
66 : CHARACTER(LEN=13) :: value
67 : CHARACTER(LEN=79) :: values
68 : CHARACTER(LEN=92) :: slice_values
69 : REAL(KIND=dp), DIMENSION(6) :: buffer, reference
70 : REAL(KIND=dp), DIMENSION(7) :: slice_buffer, slice_reference
71 :
72 2 : WRITE (value, cube_value_format) 3.33004E-101_dp
73 2 : IF (INDEX(value, "E-101") <= 0) THEN
74 0 : ERROR STOP "Cube value format must preserve three-digit negative exponents."
75 : END IF
76 :
77 : reference = [0.0_dp, -0.91458E-35_dp, 3.33004E-101_dp, &
78 2 : -1.32487E-103_dp, 1.0_dp, -1.0_dp]
79 2 : WRITE (values(1:78), cube_values_format) reference
80 2 : values(79:79) = NEW_LINE("C")
81 2 : IF (INDEX(values, "E-101") <= 0) THEN
82 0 : ERROR STOP "Cube line format must preserve E-101."
83 : END IF
84 2 : IF (INDEX(values, "E-103") <= 0) THEN
85 0 : ERROR STOP "Cube line format must preserve E-103."
86 : END IF
87 2 : CALL cube_read_values(values, buffer)
88 14 : IF (MAXVAL(ABS(buffer - reference)) > 1.0E-12_dp) THEN
89 0 : ERROR STOP "Cube reader must parse adjacent explicit-exponent values."
90 : END IF
91 :
92 14 : slice_reference(1:6) = reference
93 2 : slice_reference(7) = 0.27183E-123_dp
94 2 : WRITE (slice_values(1:78), cube_values_format) slice_reference(1:6)
95 2 : slice_values(79:79) = NEW_LINE("C")
96 2 : WRITE (slice_values(80:92), cube_value_format) slice_reference(7)
97 2 : CALL cube_read_values(slice_values, slice_buffer)
98 16 : IF (MAXVAL(ABS(slice_buffer - slice_reference)) > 1.0E-12_dp) THEN
99 0 : ERROR STOP "Cube reader must parse multi-line z-slices."
100 : END IF
101 2 : END SUBROUTINE check_cube_values
102 :
103 : ! **************************************************************************************************
104 : !> \brief Compare identical, permuted, redistributed, partial and empty PW grids.
105 : !> \param mode Layout variant; mode two changes only one rank's ordering.
106 : !> \author Thomas D. Kuehne (tkuehne@cp2k.org)
107 : ! **************************************************************************************************
108 18 : SUBROUTINE check_copy(mode)
109 : INTEGER, INTENT(IN) :: mode
110 :
111 : COMPLEX(dp), PARAMETER :: untouched = (-999.0_dp, 17.0_dp)
112 :
113 : COMPLEX(dp) :: expected
114 : INTEGER :: dimensions(1), i, j, nsource, ntarget, &
115 : owner, rank
116 : TYPE(pw_c1d_gs_type) :: destination, source
117 1854 : TYPE(pw_grid_type), TARGET :: source_grid, target_grid
118 :
119 18 : rank = para_env%mepos
120 18 : nsource = 7
121 18 : IF (mode == 6 .OR. (mode >= 7 .AND. rank == 0)) nsource = 0
122 18 : owner = rank
123 18 : IF (mode == 3 .OR. mode == 8) owner = MOD(rank + 1, para_env%num_pe)
124 18 : ntarget = nsource
125 18 : IF (mode == 4) ntarget = nsource - 1
126 18 : IF (mode == 5) ntarget = nsource + 1
127 18 : IF (mode == 8) ntarget = MERGE(0, 7, owner == 0)
128 64 : ALLOCATE (source_grid%g_hat(3, nsource), target_grid%g_hat(3, ntarget))
129 64 : ALLOCATE (source%array(nsource), destination%array(ntarget))
130 18 : source%pw_grid => source_grid
131 18 : destination%pw_grid => target_grid
132 36 : dimensions = [para_env%num_pe]
133 18 : CALL source_grid%para%group%create(para_env, 1, dimensions)
134 18 : target_grid%para%group = source_grid%para%group
135 18 : source_grid%ngpts_cut_local = nsource
136 18 : target_grid%ngpts_cut_local = ntarget
137 116 : DO i = 1, nsource
138 392 : source_grid%g_hat(:, i) = [i, rank, 0]
139 116 : source%array(i) = CMPLX(i + 11*rank, -i - 3*rank, dp)
140 : END DO
141 116 : DO i = 1, ntarget
142 98 : j = i
143 98 : IF (mode == 1 .OR. (mode == 2 .AND. rank == 0)) j = ntarget + 1 - i
144 410 : target_grid%g_hat(:, i) = [j, owner, 0]
145 : END DO
146 116 : destination%array(:) = untouched
147 18 : CALL pw_copy_match(source, destination)
148 116 : DO i = 1, ntarget
149 98 : j = target_grid%g_hat(1, i)
150 98 : expected = CMPLX(j + 11*owner, -j - 3*owner, dp)
151 98 : IF (mode == 5 .AND. j > nsource) expected = untouched
152 116 : IF (ABS(destination%array(i) - expected) > 1.e-14_dp) ERROR STOP 'Incorrect PW copy'
153 : END DO
154 18 : DEALLOCATE (source%array, destination%array)
155 18 : CALL source_grid%para%group%free()
156 36 : END SUBROUTINE check_copy
157 :
158 : ! **************************************************************************************************
159 : !> \brief Check PW multiplicity and volume normalization, and auxiliary-density Poisson routes.
160 : ! **************************************************************************************************
161 2 : SUBROUTINE check_pw_integral_and_poisson()
162 : REAL(dp), PARAMETER :: step = 5.e-4_dp, &
163 : strain_steps(4) = [step, -step, 2.0_dp*step, -2.0_dp*step]
164 :
165 : COMPLEX(dp) :: a_mode, b_mode
166 : INTEGER :: axis, ig, sign_step, span, sumtype, test
167 : REAL(dp) :: cell(3, 3), derivative, energy, energy_dv, energy_v, energy_v_dv, expected, &
168 : expected_real, kernel, reference, shifted_g(3), strain_energy(4), stress(3, 3), &
169 : stress_dv(3, 3), stress_only(3, 3), stress_v(3, 3), stress_v_dv(3, 3)
170 : TYPE(mp_comm_type) :: comm
171 14 : TYPE(pw_c1d_gs_type) :: a, b, gradient(3), gradient_v_dv(3), &
172 : potential, potential_v_dv
173 : TYPE(pw_grid_type), POINTER :: grid
174 : TYPE(pw_poisson_type) :: poisson
175 : TYPE(pw_r1d_gs_type) :: a_real, b_real
176 :
177 2 : CALL comm%set_handle(mp_comm_world%get_handle())
178 : ! Nonunit volume and an appreciable truncation factor exp(-pi).
179 2 : cell = 0.0_dp
180 8 : DO axis = 1, 3
181 8 : cell(axis, axis) = 7.0_dp**(1.0_dp/3.0_dp)
182 : END DO
183 2 : a_mode = CMPLX(0.4_dp, -0.2_dp, dp)
184 2 : b_mode = CMPLX(-0.7_dp, 0.9_dp, dp)
185 2 : NULLIFY (grid)
186 6 : DO span = FULLSPACE, HALFSPACE, -1
187 : ! Reciprocal inputs and outputs need no FFT. The production grid still
188 : ! supplies the distributed G0 ownership and half/full-space mappings.
189 4 : CALL pw_grid_create(grid, comm, cell, grid_span=span, npts=[5, 5, 5], odd=.TRUE., fft_usage=.FALSE.)
190 4 : CALL a%create(grid)
191 4 : CALL b%create(grid)
192 4 : CALL a_real%create(grid)
193 4 : CALL b_real%create(grid)
194 4 : CALL potential%create(grid)
195 4 : CALL potential_v_dv%create(grid)
196 16 : DO axis = 1, 3
197 12 : CALL gradient(axis)%create(grid)
198 16 : CALL gradient_v_dv(axis)%create(grid)
199 : END DO
200 8 : ALLOCATE (poisson%pw_pools(1))
201 4 : CALL pw_pool_create(poisson%pw_pools(1)%pool, grid)
202 4 : poisson%pw_level = 1
203 4 : poisson%parameters%solver = pw_poisson_mt
204 16 : poisson%parameters%periodic = [1, 1, 0]
205 4 : poisson%parameters%mt_alpha = 1.0_dp
206 4 : CALL pw_poisson_set(poisson, cell_hmat=cell, use_level=1)
207 12 : DO test = 1, 2
208 : ! Constant G0 only and G0 plus a Hermitian mode pair.
209 384 : a%array = 0.0_dp
210 384 : b%array = 0.0_dp
211 : expected = 0.0_dp
212 : expected_real = 0.0_dp
213 8 : reference = 0.0_dp
214 8 : IF (grid%have_g0) THEN
215 4 : a%array(1) = 2.0_dp
216 4 : b%array(1) = 3.0_dp
217 : END IF
218 8 : expected = 6.0_dp
219 8 : expected_real = 6.0_dp
220 8 : reference = 6.0_dp*mt2d_kernel_reference([0.0_dp, 0.0_dp, 0.0_dp], 3, cell(3, 3))
221 8 : IF (test /= 1) THEN
222 192 : DO ig = 1, grid%ngpts_cut_local
223 234 : IF (ANY(grid%g_hat(2:3, ig) /= 0)) CYCLE
224 12 : IF (grid%g_hat(1, ig) == 1) THEN
225 1 : a%array(ig) = a_mode
226 1 : b%array(ig) = b_mode
227 7 : ELSE IF (grid%g_hat(1, ig) == -1) THEN
228 2 : a%array(ig) = CONJG(a_mode)
229 2 : b%array(ig) = CONJG(b_mode)
230 : END IF
231 : END DO
232 4 : expected = expected + 2.0_dp*REAL(CONJG(a_mode)*b_mode, dp)
233 4 : expected_real = expected_real + 2.0_dp*REAL(a_mode, dp)*REAL(b_mode, dp)
234 16 : kernel = mt2d_kernel_reference([2.0_dp*pi/cell(1, 1), 0.0_dp, 0.0_dp], 3, cell(3, 3))
235 4 : reference = reference + 2.0_dp*REAL(CONJG(a_mode)*b_mode, dp)*kernel
236 : END IF
237 384 : a_real%array = REAL(a%array, dp)
238 384 : b_real%array = REAL(b%array, dp)
239 : ! G0 must have weight one, the mode pair weight two, and both must
240 : ! receive V=7 for an integral. just_sum retains the unscaled oracle.
241 24 : DO sumtype = do_accurate_sum, do_standard_sum
242 16 : CALL check_error(pw_integral_ab(a, b, sumtype=sumtype), 7.0_dp*expected, 3.e-14_dp)
243 16 : CALL check_error(pw_integral_ab(a, b, sumtype=sumtype, just_sum=.TRUE.), expected, 3.e-14_dp)
244 16 : CALL check_error(pw_integral_ab(a_real, b_real, sumtype=sumtype), 7.0_dp*expected_real, 3.e-14_dp)
245 : CALL check_error(pw_integral_ab(a_real, b_real, sumtype=sumtype, just_sum=.TRUE.), &
246 16 : expected_real, 3.e-14_dp)
247 16 : CALL check_error(pw_integral_ab(a, b_real, sumtype=sumtype), 7.0_dp*expected_real, 3.e-14_dp)
248 24 : CALL check_error(pw_integral_ab(a_real, b, sumtype=sumtype), 7.0_dp*expected_real, 3.e-14_dp)
249 : END DO
250 8 : CALL pw_poisson_solve(poisson, a, energy, h_stress=stress, aux_density=b)
251 8 : CALL pw_poisson_solve(poisson, a, energy_dv, dvhartree=gradient, h_stress=stress_dv, aux_density=b)
252 8 : CALL pw_poisson_solve(poisson, a, energy_v, potential, h_stress=stress_v, aux_density=b)
253 : CALL pw_poisson_solve(poisson, a, energy_v_dv, potential_v_dv, dvhartree=gradient_v_dv, &
254 8 : h_stress=stress_v_dv, aux_density=b)
255 8 : CALL pw_poisson_solve(poisson, a, h_stress=stress_only, aux_density=b)
256 8 : CALL check_error(energy, 0.5_dp*7.0_dp*reference, 3.e-13_dp)
257 8 : CALL check_error(energy_dv, energy, 3.e-13_dp)
258 8 : CALL check_error(energy_v, energy, 3.e-13_dp)
259 8 : CALL check_error(energy_v_dv, energy, 3.e-13_dp)
260 104 : CALL check_error(MAXVAL(ABS(stress_dv - stress)), 0.0_dp, 3.e-11_dp)
261 104 : CALL check_error(MAXVAL(ABS(stress_v - stress)), 0.0_dp, 3.e-11_dp)
262 104 : CALL check_error(MAXVAL(ABS(stress_v_dv - stress)), 0.0_dp, 3.e-11_dp)
263 104 : CALL check_error(MAXVAL(ABS(stress_only - stress)), 0.0_dp, 3.e-11_dp)
264 8 : IF (test == 2) THEN
265 : ! Fixed Fourier charges Q=V*rho give E(epsilon)=V/2/(1+epsilon)
266 : ! times the strained kernel sum. This absolute production-virial
267 : ! oracle includes both G0/volume and the nonzero-mode kernel term.
268 20 : DO sign_step = 1, 4
269 64 : shifted_g = [2.0_dp*pi/cell(1, 1)/(1.0_dp + strain_steps(sign_step)), 0.0_dp, 0.0_dp]
270 : reference = 6.0_dp*mt2d_kernel_reference([0.0_dp, 0.0_dp, 0.0_dp], 3, cell(3, 3)) &
271 16 : + 2.0_dp*REAL(CONJG(a_mode)*b_mode, dp)*mt2d_kernel_reference(shifted_g, 3, cell(3, 3))
272 20 : strain_energy(sign_step) = 0.5_dp*7.0_dp*reference/(1.0_dp + strain_steps(sign_step))
273 : END DO
274 : derivative = (8.0_dp*(strain_energy(1) - strain_energy(2)) &
275 4 : - strain_energy(3) + strain_energy(4))/(12.0_dp*step)
276 4 : CALL check_error(stress(1, 1) + derivative, 0.0_dp, 3.e-9_dp)
277 : END IF
278 404 : DO ig = 1, grid%ngpts_cut_local
279 376 : kernel = mt2d_kernel_reference(grid%g(:, ig), 3, cell(3, 3))
280 376 : CALL check_error(ABS(potential%array(ig) - kernel*a%array(ig)), 0.0_dp, 3.e-13_dp)
281 376 : CALL check_error(ABS(potential_v_dv%array(ig) - kernel*a%array(ig)), 0.0_dp, 3.e-13_dp)
282 1512 : DO axis = 1, 3
283 : CALL check_error(ABS(gradient(axis)%array(ig) &
284 1128 : - CMPLX(0.0_dp, grid%g(axis, ig), dp)*kernel*a%array(ig)), 0.0_dp, 3.e-13_dp)
285 : CALL check_error(ABS(gradient_v_dv(axis)%array(ig) &
286 1504 : - CMPLX(0.0_dp, grid%g(axis, ig), dp)*kernel*a%array(ig)), 0.0_dp, 3.e-13_dp)
287 : END DO
288 : END DO
289 : END DO
290 4 : CALL poisson%release()
291 4 : CALL a%release()
292 4 : CALL b%release()
293 4 : CALL a_real%release()
294 4 : CALL b_real%release()
295 4 : CALL potential%release()
296 4 : CALL potential_v_dv%release()
297 16 : DO axis = 1, 3
298 12 : CALL gradient(axis)%release()
299 16 : CALL gradient_v_dv(axis)%release()
300 : END DO
301 6 : CALL pw_grid_release(grid)
302 : END DO
303 34 : END SUBROUTINE check_pw_integral_and_poisson
304 :
305 : ! **************************************************************************************************
306 : !> \brief Compare the MT2D kernel and in-plane virial with independent strained-cell energies.
307 : ! **************************************************************************************************
308 2 : SUBROUTINE check_mt2d_kernel()
309 : REAL(dp), PARAMETER :: step = 3.e-4_dp, shifts(4) = [step, -step, 2.0_dp*step, -2.0_dp*step]
310 :
311 : INTEGER :: alpha, beta, ip, iz, jp, miller(3), &
312 : mode, modes(3, 5), sign_step
313 : REAL(dp) :: derivative, g(3), h_inv(3, 3), kernel, &
314 : lengths(3), shifted_g(3), &
315 : shifted_kernel(4), strain(3, 3)
316 :
317 : ! Rotate the Cartesian vacuum axis. A shear within the periodic plane is
318 : ! permitted; only normal/periodic mixing is excluded by the MT2D kernel.
319 8 : DO iz = 1, 3
320 6 : ip = MOD(iz, 3) + 1
321 6 : jp = MOD(iz + 1, 3) + 1
322 6 : lengths = [9.0_dp, 13.0_dp, 17.0_dp]
323 6 : h_inv = 0.0_dp
324 24 : DO alpha = 1, 3
325 24 : h_inv(alpha, alpha) = 1.0_dp/lengths(alpha)
326 : END DO
327 6 : h_inv(ip, jp) = -1.4_dp/(lengths(ip)*lengths(jp))
328 6 : modes = 0
329 : ! G=0, both odd and even normal-only modes, and several modes with
330 : ! appreciable exponential truncation reject the periodic K**2 formula.
331 18 : modes(iz, 2:3) = [1, 2]
332 6 : modes(ip, 4) = 1
333 24 : modes(:, 5) = [1, 2, 1]
334 38 : DO mode = 1, SIZE(modes, 2)
335 120 : miller = modes(:, mode)
336 240 : g = 2.0_dp*pi*MATMUL(TRANSPOSE(h_inv), REAL(miller, dp))
337 30 : CALL mt2d_kernel_strain(g, iz, lengths(iz), kernel, strain)
338 30 : CALL check_error(kernel, mt2d_kernel_reference(g, iz, lengths(iz)), 3.e-14_dp)
339 120 : CALL check_error(MAXVAL(ABS(strain(iz, :))), 0.0_dp, 0.0_dp)
340 120 : CALL check_error(MAXVAL(ABS(strain(:, iz))), 0.0_dp, 0.0_dp)
341 246 : IF (mode <= 3) CALL check_error(MAXVAL(ABS(strain)), 0.0_dp, 0.0_dp)
342 156 : DO alpha = 1, 3
343 90 : IF (alpha == iz) CYCLE
344 270 : DO beta = 1, 3
345 180 : IF (beta == iz) CYCLE
346 600 : DO sign_step = 1, 4
347 : ! g'=(I+epsilon_ab)^(-T) g, evaluated independently of dK.
348 480 : shifted_g = g
349 480 : IF (alpha == beta) THEN
350 240 : shifted_g(beta) = g(beta)/(1.0_dp + shifts(sign_step))
351 : ELSE
352 240 : shifted_g(beta) = g(beta) - shifts(sign_step)*g(alpha)
353 : END IF
354 600 : shifted_kernel(sign_step) = mt2d_kernel_reference(shifted_g, iz, lengths(iz))
355 : END DO
356 : derivative = (8.0_dp*(shifted_kernel(1) - shifted_kernel(2)) &
357 120 : - shifted_kernel(3) + shifted_kernel(4))/(12.0_dp*step)
358 270 : CALL check_error(strain(alpha, beta), derivative, 3.e-9_dp)
359 : END DO
360 : END DO
361 : END DO
362 : END DO
363 2 : END SUBROUTINE check_mt2d_kernel
364 :
365 : ! **************************************************************************************************
366 : !> \brief Evaluate the original MT2D influence kernel independently of its strain derivative.
367 : !> \param g ...
368 : !> \param iz ...
369 : !> \param slab_size ...
370 : !> \return ...
371 : ! **************************************************************************************************
372 930 : FUNCTION mt2d_kernel_reference(g, iz, slab_size) RESULT(kernel)
373 : REAL(dp), INTENT(IN) :: g(3)
374 : INTEGER, INTENT(IN) :: iz
375 : REAL(dp), INTENT(IN) :: slab_size
376 : REAL(dp) :: kernel
377 :
378 : REAL(dp) :: g2, g_parallel(3), radius
379 :
380 3720 : g2 = SUM(g**2)
381 930 : IF (g2 == 0.0_dp) THEN
382 130 : kernel = pi*slab_size**2/2.0_dp
383 : ELSE
384 800 : g_parallel = g
385 800 : g_parallel(iz) = 0.0_dp
386 800 : radius = slab_size/2.0_dp
387 3200 : kernel = fourpi/g2 - fourpi/g2*COS(g(iz)*radius)*EXP(-NORM2(g_parallel)*radius)
388 : END IF
389 930 : END FUNCTION mt2d_kernel_reference
390 :
391 : ! **************************************************************************************************
392 : !> \brief Require the scaled error to be finite and within the supplied tolerance.
393 : !> \param value ...
394 : !> \param reference ...
395 : !> \param tolerance ...
396 : ! **************************************************************************************************
397 3400 : SUBROUTINE check_error(value, reference, tolerance)
398 : REAL(dp), INTENT(IN) :: value, reference, tolerance
399 :
400 : REAL(dp) :: error
401 :
402 3400 : error = ABS(value - reference)/MAX(1.0_dp, ABS(reference))
403 3400 : max_error = MAX(max_error, error)
404 3400 : IF (.NOT. (error <= tolerance)) THEN
405 0 : PRINT *, 'Coulomb check: value, reference, error, tolerance:', value, reference, error, tolerance
406 0 : ERROR STOP 'PW integral or MT2D kernel check failed'
407 : END IF
408 3400 : END SUBROUTINE check_error
409 :
410 : END PROGRAM pw_unittest
|
