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 3970 : PROGRAM kpoint_lattice_fft_unittest
8 2 : USE fft_tools, ONLY: finalize_fft,&
9 : init_fft
10 : USE kinds, ONLY: dp
11 : USE kpoint_lattice_fft, ONLY: cell_to_k_grid_fft,&
12 : regular_kpoint_grid
13 : USE mathconstants, ONLY: gaussi,&
14 : twopi,&
15 : z_zero
16 : USE message_passing, ONLY: mp_comm_type,&
17 : mp_world_finalize,&
18 : mp_world_init
19 :
20 : IMPLICIT NONE
21 :
22 : INTEGER, PARAMETER :: ncell = 31, nkp = 64, nselected = 3
23 : COMPLEX(KIND=dp), DIMENSION(2, 3, nkp) :: direct, transformed
24 : COMPLEX(KIND=dp), DIMENSION(2, 3, nselected) :: selected_transformed
25 : INTEGER :: d, icell, ik, ix, iy, iz
26 : INTEGER, DIMENSION(3) :: nkp_grid
27 : INTEGER, DIMENSION(3, ncell) :: index_to_cell
28 : INTEGER, DIMENSION(nselected) :: selected_kpoints
29 : LOGICAL :: used_fft
30 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
31 : REAL(KIND=dp), DIMENSION(3, nkp) :: xkp, xkp_bad
32 : REAL(KIND=dp), DIMENSION(2, 3, ncell) :: values_rs
33 : TYPE(mp_comm_type) :: mp_comm
34 :
35 2 : CALL mp_world_init(mp_comm)
36 : CALL init_fft("FFTSG", alltoall=.FALSE., fftsg_sizes=.TRUE., pool_limit=4, &
37 2 : wisdom_file="", plan_style=1)
38 :
39 2 : nkp_grid = [4, 4, 4]
40 : hmat = RESHAPE([4.0_dp, 0.1_dp, 0.2_dp, &
41 : 0.0_dp, 5.0_dp, 0.3_dp, &
42 2 : 0.0_dp, 0.0_dp, 6.0_dp], [3, 3])
43 :
44 2 : ik = 0
45 10 : DO iz = 0, nkp_grid(3) - 1
46 42 : DO iy = 0, nkp_grid(2) - 1
47 168 : DO ix = 0, nkp_grid(1) - 1
48 128 : ik = ik + 1
49 : xkp(:, ik) = [(REAL(ix, dp) + 0.5_dp)/REAL(nkp_grid(1), dp), &
50 : REAL(iy, dp)/REAL(nkp_grid(2), dp), &
51 544 : (REAL(iz, dp) + 0.25_dp)/REAL(nkp_grid(3), dp)]
52 : END DO
53 : END DO
54 : END DO
55 2 : IF (.NOT. regular_kpoint_grid(xkp, nkp_grid)) ERROR STOP "Regular k-point grid not recognized"
56 :
57 64 : DO icell = 1, ncell
58 : index_to_cell(:, icell) = [MODULO(3*icell + 1, 11) - 5, &
59 : MODULO(5*icell + 2, 9) - 4, &
60 248 : MODULO(7*icell + 3, 7) - 3]
61 250 : DO d = 1, 3
62 186 : values_rs(1, d, icell) = SIN(0.13_dp*REAL(icell + 2*d, dp))
63 248 : values_rs(2, d, icell) = COS(0.17_dp*REAL(2*icell - d, dp))
64 : END DO
65 : END DO
66 :
67 2 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, transformed, used_fft)
68 2 : IF (.NOT. used_fft) ERROR STOP "Complete regular grid did not use lattice FFT"
69 2 : CALL direct_reference(values_rs, index_to_cell, xkp, direct)
70 1282 : IF (MAXVAL(ABS(transformed - direct)) > 2.0E-12_dp) THEN
71 0 : ERROR STOP "Shifted lattice FFT differs from direct transform"
72 : END IF
73 :
74 2 : selected_kpoints = [2, 31, 64]
75 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, selected_transformed, &
76 2 : used_fft, selected_kpoints=selected_kpoints)
77 2 : IF (.NOT. used_fft) ERROR STOP "Selected distributed output did not use lattice FFT"
78 62 : IF (MAXVAL(ABS(selected_transformed - direct(:, :, selected_kpoints))) > 2.0E-12_dp) THEN
79 0 : ERROR STOP "Selected lattice-FFT output differs from the full transform"
80 : END IF
81 :
82 2 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, transformed, used_fft, 2, hmat)
83 2 : IF (.NOT. used_fft) ERROR STOP "Derivative did not use lattice FFT"
84 2 : CALL direct_reference(values_rs, index_to_cell, xkp, direct, 2, hmat)
85 1282 : IF (MAXVAL(ABS(transformed - direct)) > 5.0E-12_dp) THEN
86 0 : ERROR STOP "Lattice-FFT derivative differs from direct transform"
87 : END IF
88 :
89 2 : nkp_grid = [4, 4, 1]
90 2 : ik = 0
91 10 : DO iy = 0, nkp_grid(2) - 1
92 42 : DO ix = 0, nkp_grid(1) - 1
93 32 : ik = ik + 1
94 : xkp(:, ik) = [(REAL(ix, dp) + 0.5_dp)/REAL(nkp_grid(1), dp), &
95 136 : REAL(iy, dp)/REAL(nkp_grid(2), dp), 0.0_dp]
96 : END DO
97 : END DO
98 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp(:, 1:ik), nkp_grid, &
99 2 : transformed(:, :, 1:ik), used_fft)
100 2 : IF (.NOT. used_fft) ERROR STOP "Singleton grid direction did not use padded lattice FFT"
101 2 : CALL direct_reference(values_rs, index_to_cell, xkp(:, 1:ik), direct(:, :, 1:ik))
102 322 : IF (MAXVAL(ABS(transformed(:, :, 1:ik) - direct(:, :, 1:ik))) > 2.0E-12_dp) THEN
103 0 : ERROR STOP "Padded lattice FFT differs from direct transform"
104 : END IF
105 :
106 2 : nkp_grid = [4, 4, 4]
107 :
108 2 : xkp_bad = xkp
109 2 : xkp_bad(1, nkp) = xkp_bad(1, nkp) + 0.03125_dp
110 2 : CALL cell_to_k_grid_fft(values_rs, index_to_cell, xkp_bad, nkp_grid, transformed, used_fft)
111 2 : IF (used_fft) ERROR STOP "Irregular k-point list incorrectly used lattice FFT"
112 2 : CALL direct_reference(values_rs, index_to_cell, xkp_bad, direct)
113 1282 : IF (MAXVAL(ABS(transformed - direct)) > 1.0E-14_dp) THEN
114 0 : ERROR STOP "Direct lattice-transform fallback is inconsistent"
115 : END IF
116 :
117 2 : CALL finalize_fft(mp_comm, "")
118 2 : CALL mp_world_finalize()
119 :
120 : CONTAINS
121 :
122 : ! **************************************************************************************************
123 : !> \brief ...
124 : !> \param values_rs ...
125 : !> \param index_to_cell ...
126 : !> \param xkp ...
127 : !> \param values_k ...
128 : !> \param deriv_direction ...
129 : !> \param hmat ...
130 : ! **************************************************************************************************
131 8 : SUBROUTINE direct_reference(values_rs, index_to_cell, xkp, values_k, deriv_direction, hmat)
132 :
133 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: values_rs
134 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
135 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
136 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_k
137 : INTEGER, INTENT(IN), OPTIONAL :: deriv_direction
138 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
139 : OPTIONAL :: hmat
140 :
141 : COMPLEX(KIND=dp) :: factor
142 : INTEGER :: icell, ik
143 : REAL(KIND=dp), DIMENSION(3) :: cell_vector
144 :
145 4168 : values_k = z_zero
146 424 : DO ik = 1, SIZE(xkp, 2)
147 13320 : DO icell = 1, SIZE(values_rs, 3)
148 51584 : factor = EXP(gaussi*twopi*SUM(xkp(:, ik)*REAL(index_to_cell(:, icell), KIND=dp)))
149 12896 : IF (PRESENT(deriv_direction)) THEN
150 63488 : cell_vector = MATMUL(hmat, REAL(index_to_cell(:, icell), KIND=dp))
151 3968 : factor = factor*gaussi*cell_vector(deriv_direction)
152 : END IF
153 129376 : values_k(:, :, ik) = values_k(:, :, ik) + factor*values_rs(:, :, icell)
154 : END DO
155 : END DO
156 :
157 8 : END SUBROUTINE direct_reference
158 :
159 : END PROGRAM kpoint_lattice_fft_unittest
|
