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 : ! **************************************************************************************************
9 : !> \brief Batched lattice Fourier transforms for complete regular k-point grids.
10 : !>
11 : !> Real-space image cells are folded modulo the reciprocal grid. Shifted Monkhorst-Pack
12 : !> meshes are handled by a phase twist before the FFT. Non-regular or incomplete k-point
13 : !> lists transparently use the direct phase sum.
14 : ! **************************************************************************************************
15 : MODULE kpoint_lattice_fft
16 : USE fft_tools, ONLY: BWFFT,&
17 : FFT_RADIX_NEXT,&
18 : fft3d,&
19 : fft_alloc,&
20 : fft_dealloc,&
21 : fft_radix_operations,&
22 : fft_type
23 : USE kinds, ONLY: dp
24 : USE mathconstants, ONLY: gaussi,&
25 : twopi,&
26 : z_zero
27 : #include "./base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 :
31 : PRIVATE
32 :
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_lattice_fft'
34 :
35 : PUBLIC :: cell_to_k_grid_fft, &
36 : regular_kpoint_grid
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief Test and map a complete, uniformly shifted reciprocal grid.
42 : !> \param xkp reciprocal coordinates, in units of reciprocal lattice vectors
43 : !> \param nkp_grid regular-grid dimensions
44 : !> \param grid_index FFT-grid index for every k point
45 : !> \param k_offset common reciprocal-coordinate offset
46 : !> \return true for a complete regular grid
47 : ! **************************************************************************************************
48 28 : LOGICAL FUNCTION regular_kpoint_grid(xkp, nkp_grid, grid_index, k_offset) RESULT(regular)
49 :
50 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
51 : INTEGER, DIMENSION(3), INTENT(IN) :: nkp_grid
52 : INTEGER, DIMENSION(:, :), INTENT(OUT), OPTIONAL :: grid_index
53 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: k_offset
54 :
55 : REAL(KIND=dp), PARAMETER :: map_tolerance = 2.0E-10_dp
56 :
57 : INTEGER :: d, ik, m, nkp
58 : INTEGER, DIMENSION(3) :: index
59 28 : LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: occupied
60 : REAL(KIND=dp) :: offset(3), scaled
61 :
62 28 : regular = .FALSE.
63 112 : IF (SIZE(xkp, 1) < 3 .OR. ANY(nkp_grid <= 0)) RETURN
64 28 : nkp = SIZE(xkp, 2)
65 112 : IF (PRODUCT(nkp_grid) /= nkp) RETURN
66 28 : IF (PRESENT(grid_index)) THEN
67 18 : IF (SIZE(grid_index, 1) < 3 .OR. SIZE(grid_index, 2) < nkp) RETURN
68 : END IF
69 :
70 112 : DO d = 1, 3
71 84 : scaled = MODULO(REAL(nkp_grid(d), KIND=dp)*xkp(d, 1), 1.0_dp)
72 84 : IF (ABS(scaled) < map_tolerance .OR. ABS(scaled - 1.0_dp) < map_tolerance) scaled = 0.0_dp
73 112 : offset(d) = scaled/REAL(nkp_grid(d), KIND=dp)
74 : END DO
75 :
76 140 : ALLOCATE (occupied(nkp_grid(1), nkp_grid(2), nkp_grid(3)), source=.FALSE.)
77 1628 : DO ik = 1, nkp
78 6406 : DO d = 1, 3
79 4806 : scaled = REAL(nkp_grid(d), KIND=dp)*(xkp(d, ik) - offset(d))
80 4806 : m = NINT(scaled)
81 4806 : IF (ABS(scaled - REAL(m, KIND=dp)) > map_tolerance) THEN
82 2 : DEALLOCATE (occupied)
83 2 : RETURN
84 : END IF
85 6404 : INDEX(d) = MODULO(m, nkp_grid(d)) + 1
86 : END DO
87 1600 : IF (occupied(INDEX(1), INDEX(2), INDEX(3))) THEN
88 0 : DEALLOCATE (occupied)
89 0 : RETURN
90 : END IF
91 1600 : occupied(INDEX(1), INDEX(2), INDEX(3)) = .TRUE.
92 4506 : IF (PRESENT(grid_index)) grid_index(1:3, ik) = index
93 : END DO
94 2084 : regular = ALL(occupied)
95 26 : DEALLOCATE (occupied)
96 :
97 26 : IF (regular .AND. PRESENT(k_offset)) k_offset = offset
98 :
99 : END FUNCTION regular_kpoint_grid
100 :
101 : ! **************************************************************************************************
102 : !> \brief Transform a batch of real matrices from image cells to every supplied k point.
103 : !> \param values_rs real-space matrices; the last dimension enumerates image cells
104 : !> \param index_to_cell integer lattice vector for every image cell
105 : !> \param xkp reciprocal coordinates of all requested k points
106 : !> \param nkp_grid dimensions of the candidate regular reciprocal grid
107 : !> \param values_k complex matrices at every k point
108 : !> \param used_fft reports whether the regular-grid FFT path was used
109 : !> \param deriv_direction optional Cartesian k derivative (1, 2, or 3)
110 : !> \param hmat direct-lattice cell matrix, required for a derivative
111 : !> \param selected_kpoints optional indices of k points to retain in values_k
112 : ! **************************************************************************************************
113 36 : SUBROUTINE cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, values_k, &
114 18 : used_fft, deriv_direction, hmat, selected_kpoints)
115 :
116 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: values_rs
117 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
118 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
119 : INTEGER, DIMENSION(3), INTENT(IN) :: nkp_grid
120 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_k
121 : LOGICAL, INTENT(OUT), OPTIONAL :: used_fft
122 : INTEGER, INTENT(IN), OPTIONAL :: deriv_direction
123 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
124 : OPTIONAL :: hmat
125 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: selected_kpoints
126 :
127 18 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cell_factor
128 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
129 18 : POINTER :: fft_in, fft_out
130 : INTEGER :: attempt, d, handle, i, icell, ik, &
131 : ik_out, j, nout, radix_length, stat
132 18 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: cell_index, grid_index
133 : INTEGER, DIMENSION(3) :: n, nfft
134 : LOGICAL :: compatible, regular
135 : REAL(KIND=dp) :: arg
136 : REAL(KIND=dp), DIMENSION(3) :: cell_vector, k_offset
137 :
138 18 : CALL timeset("cell_to_k_grid_fft", handle)
139 :
140 18 : IF (SIZE(index_to_cell, 1) < 3 .OR. &
141 : SIZE(index_to_cell, 2) /= SIZE(values_rs, 3)) THEN
142 0 : CPABORT("Inconsistent real-space image-cell mapping in lattice FFT")
143 : END IF
144 18 : nout = SIZE(xkp, 2)
145 18 : IF (PRESENT(selected_kpoints)) THEN
146 2 : nout = SIZE(selected_kpoints)
147 14 : IF (ANY(selected_kpoints < 1) .OR. ANY(selected_kpoints > SIZE(xkp, 2))) THEN
148 0 : CPABORT("Selected lattice-FFT k-point index is out of range")
149 : END IF
150 : END IF
151 : IF (SIZE(values_k, 1) /= SIZE(values_rs, 1) .OR. &
152 18 : SIZE(values_k, 2) /= SIZE(values_rs, 2) .OR. &
153 : SIZE(values_k, 3) /= nout) THEN
154 0 : CPABORT("Inconsistent input and output matrix batches in lattice FFT")
155 : END IF
156 18 : IF (PRESENT(deriv_direction)) THEN
157 2 : IF (.NOT. PRESENT(hmat)) CPABORT("Lattice-FFT derivative requested without a cell matrix")
158 2 : IF (deriv_direction < 1 .OR. deriv_direction > 3) THEN
159 0 : CPABORT("Lattice-FFT derivative direction must be 1, 2, or 3")
160 : END IF
161 : END IF
162 :
163 54 : ALLOCATE (grid_index(3, SIZE(xkp, 2)))
164 18 : regular = regular_kpoint_grid(xkp, nkp_grid, grid_index, k_offset)
165 18 : n = nkp_grid
166 18 : nfft = n
167 18 : IF (regular .AND. fft_type /= 3) THEN
168 32 : DO d = 1, 3
169 : ! FFTSG rejects some short or non-radix lengths. Embed the N-point
170 : ! transform in a supported integer multiple L: real-space samples
171 : ! occupy every L/N-th slot, while the first N output frequencies
172 : ! are exactly the desired N-point DFT. This keeps common 4x4x1
173 : ! and 2D/1D Monkhorst-Pack grids on the FFT path.
174 30 : compatible = .FALSE.
175 30 : DO attempt = 0, 15
176 30 : IF (nfft(d) >= 3) THEN
177 26 : CALL fft_radix_operations(nfft(d), radix_length, FFT_RADIX_NEXT)
178 26 : IF (radix_length == nfft(d)) THEN
179 : compatible = .TRUE.
180 : EXIT
181 : END IF
182 : END IF
183 6 : nfft(d) = nfft(d) + n(d)
184 : END DO
185 32 : IF (.NOT. compatible) regular = .FALSE.
186 : END DO
187 : END IF
188 16 : IF (.NOT. regular) THEN
189 : CALL direct_cell_to_k(values_rs, index_to_cell, xkp, values_k, deriv_direction, hmat, &
190 4 : selected_kpoints)
191 2 : IF (PRESENT(used_fft)) used_fft = .FALSE.
192 2 : DEALLOCATE (grid_index)
193 2 : CALL timestop(handle)
194 : RETURN
195 : END IF
196 :
197 16 : NULLIFY (fft_in, fft_out)
198 16 : CALL fft_alloc(fft_in, nfft)
199 16 : CALL fft_alloc(fft_out, nfft)
200 80 : ALLOCATE (cell_index(3, SIZE(values_rs, 3)), cell_factor(SIZE(values_rs, 3)))
201 1696 : DO icell = 1, SIZE(values_rs, 3)
202 6720 : DO d = 1, 3
203 6720 : cell_index(d, icell) = MODULO(index_to_cell(d, icell), n(d))*(nfft(d)/n(d)) + 1
204 : END DO
205 6720 : arg = SUM(k_offset*REAL(index_to_cell(1:3, icell), KIND=dp))
206 1680 : cell_factor(icell) = EXP(gaussi*twopi*arg)
207 1696 : IF (PRESENT(deriv_direction)) THEN
208 992 : cell_vector = MATMUL(hmat, REAL(index_to_cell(1:3, icell), KIND=dp))
209 62 : cell_factor(icell) = cell_factor(icell)*gaussi*cell_vector(deriv_direction)
210 : END IF
211 : END DO
212 :
213 16 : stat = 0
214 48 : DO j = 1, SIZE(values_rs, 2)
215 3680 : DO i = 1, SIZE(values_rs, 1)
216 3648 : fft_in = z_zero
217 649536 : DO icell = 1, SIZE(values_rs, 3)
218 : fft_in(cell_index(1, icell), cell_index(2, icell), cell_index(3, icell)) = &
219 : fft_in(cell_index(1, icell), cell_index(2, icell), cell_index(3, icell)) + &
220 649536 : cell_factor(icell)*values_rs(i, j, icell)
221 : END DO
222 3648 : CALL fft3d(BWFFT, nfft, fft_in, fft_out, status=stat)
223 3648 : IF (stat /= 0) EXIT
224 235844 : DO ik_out = 1, nout
225 232164 : ik = ik_out
226 232164 : IF (PRESENT(selected_kpoints)) ik = selected_kpoints(ik_out)
227 : values_k(i, j, ik_out) = &
228 235812 : fft_out(grid_index(1, ik), grid_index(2, ik), grid_index(3, ik))
229 : END DO
230 : END DO
231 48 : IF (stat /= 0) EXIT
232 : END DO
233 :
234 16 : IF (stat /= 0) THEN
235 : CALL direct_cell_to_k(values_rs, index_to_cell, xkp, values_k, deriv_direction, hmat, &
236 0 : selected_kpoints)
237 0 : regular = .FALSE.
238 : END IF
239 16 : IF (PRESENT(used_fft)) used_fft = regular
240 :
241 16 : CALL fft_dealloc(fft_in)
242 16 : CALL fft_dealloc(fft_out)
243 16 : DEALLOCATE (cell_factor, cell_index, grid_index)
244 16 : CALL timestop(handle)
245 :
246 36 : END SUBROUTINE cell_to_k_grid_fft
247 :
248 : ! **************************************************************************************************
249 : !> \brief Direct phase-sum fallback for arbitrary reciprocal-point lists.
250 : !> \param values_rs ...
251 : !> \param index_to_cell ...
252 : !> \param xkp ...
253 : !> \param values_k ...
254 : !> \param deriv_direction ...
255 : !> \param hmat ...
256 : !> \param selected_kpoints ...
257 : ! **************************************************************************************************
258 4 : SUBROUTINE direct_cell_to_k(values_rs, index_to_cell, xkp, values_k, &
259 2 : deriv_direction, hmat, selected_kpoints)
260 :
261 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: values_rs
262 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
263 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
264 : COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_k
265 : INTEGER, INTENT(IN), OPTIONAL :: deriv_direction
266 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
267 : OPTIONAL :: hmat
268 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: selected_kpoints
269 :
270 : COMPLEX(KIND=dp) :: factor
271 : INTEGER :: icell, ik, ik_out
272 : REAL(KIND=dp) :: arg
273 : REAL(KIND=dp), DIMENSION(3) :: cell_vector
274 :
275 1282 : values_k = z_zero
276 130 : DO ik_out = 1, SIZE(values_k, 3)
277 128 : ik = ik_out
278 128 : IF (PRESENT(selected_kpoints)) ik = selected_kpoints(ik_out)
279 4098 : DO icell = 1, SIZE(values_rs, 3)
280 15872 : arg = SUM(xkp(1:3, ik)*REAL(index_to_cell(1:3, icell), KIND=dp))
281 3968 : factor = EXP(gaussi*twopi*arg)
282 3968 : IF (PRESENT(deriv_direction)) THEN
283 0 : cell_vector = MATMUL(hmat, REAL(index_to_cell(1:3, icell), KIND=dp))
284 0 : factor = factor*gaussi*cell_vector(deriv_direction)
285 : END IF
286 39808 : values_k(:, :, ik_out) = values_k(:, :, ik_out) + factor*values_rs(:, :, icell)
287 : END DO
288 : END DO
289 :
290 2 : END SUBROUTINE direct_cell_to_k
291 :
292 : END MODULE kpoint_lattice_fft
|
