(git:6344619)
Loading...
Searching...
No Matches
kpoint_lattice_fft.F
Go to the documentation of this file.
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 on regular k-point grids.
10!>
11!> Image cells are folded modulo the reciprocal grid; shifted meshes use a phase twist.
12!> Forward transforms require a complete grid. The inverse transform can also embed
13!> an explicitly allowed subset with its supplied weights and zero missing points.
14!> Unsupported grids retain the direct phase-sum fallback.
15! **************************************************************************************************
16MODULE kpoint_lattice_fft
17 USE fft_lib, ONLY: fft_supports_arbitrary_lengths
18 USE fft_tools, ONLY: bwfft,&
19 fft_radix_next,&
20 fwfft,&
21 fft3d,&
22 fft_alloc,&
23 fft_dealloc,&
24 fft_radix_operations
25 USE kinds, ONLY: dp
26 USE mathconstants, ONLY: gaussi,&
27 twopi,&
28 z_zero
29#include "./base/base_uses.f90"
30
31 IMPLICIT NONE
32
33 PRIVATE
34
35 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_lattice_fft'
36
37 PUBLIC :: cell_to_k_grid_fft, &
38 k_grid_to_cell_fft, &
39 k_grid_to_cell_prepare, k_grid_to_cell_execute, k_grid_to_cell_release, &
40 k_grid_to_cell_work_type, &
41 lattice_fft_shape, &
42 regular_kpoint_grid
43
44 ! Call-scoped mathematical maps and buffers. Backend plans remain owned by fft_tools.
45 TYPE k_grid_to_cell_work_type
46 PRIVATE
47 LOGICAL :: regular = .false.
48 INTEGER :: nfft(3) = 0
49 INTEGER, ALLOCATABLE :: cells(:, :), grid_index(:, :), cell_index(:, :)
50 REAL(kind=dp), ALLOCATABLE :: xkp(:, :), weight(:)
51 COMPLEX(KIND=dp), ALLOCATABLE :: cell_phase(:)
52 COMPLEX(KIND=dp), CONTIGUOUS, POINTER :: fft_in(:, :, :) => null(), fft_out(:, :, :) => null()
53 END TYPE k_grid_to_cell_work_type
54
55CONTAINS
56
57! **************************************************************************************************
58!> \brief Test and map a uniformly shifted reciprocal grid or an explicitly allowed subset.
59!> \param xkp reciprocal coordinates, in units of reciprocal lattice vectors
60!> \param nkp_grid regular-grid dimensions
61!> \param grid_index FFT-grid index for every k point
62!> \param k_offset common reciprocal-coordinate offset
63!> \param allow_incomplete accept a subset of a regular grid (missing points are zero)
64!> \return true for a regular grid, complete unless allow_incomplete is set
65! **************************************************************************************************
66 LOGICAL FUNCTION regular_kpoint_grid(xkp, nkp_grid, grid_index, k_offset, allow_incomplete) RESULT(regular)
67
68 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
69 INTEGER, DIMENSION(3), INTENT(IN) :: nkp_grid
70 INTEGER, DIMENSION(:, :), INTENT(OUT), OPTIONAL :: grid_index
71 REAL(kind=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: k_offset
72 LOGICAL, INTENT(IN), OPTIONAL :: allow_incomplete
73
74 REAL(kind=dp), PARAMETER :: map_tolerance = 2.0e-10_dp
75
76 INTEGER :: d, ik, m, nkp
77 INTEGER, DIMENSION(3) :: index
78 LOGICAL :: require_complete
79 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: occupied
80 REAL(kind=dp) :: offset(3), scaled
81
82 regular = .false.
83 IF (SIZE(xkp, 1) < 3 .OR. any(nkp_grid <= 0)) RETURN
84 nkp = SIZE(xkp, 2)
85 require_complete = .true.
86 IF (PRESENT(allow_incomplete)) require_complete = .NOT. allow_incomplete
87 IF (nkp < 1 .OR. nkp > product(nkp_grid)) RETURN
88 IF (require_complete .AND. product(nkp_grid) /= nkp) RETURN
89 IF (PRESENT(grid_index)) THEN
90 IF (SIZE(grid_index, 1) < 3 .OR. SIZE(grid_index, 2) < nkp) RETURN
91 END IF
92
93 DO d = 1, 3
94 scaled = modulo(real(nkp_grid(d), kind=dp)*xkp(d, 1), 1.0_dp)
95 IF (abs(scaled) < map_tolerance .OR. abs(scaled - 1.0_dp) < map_tolerance) scaled = 0.0_dp
96 offset(d) = scaled/real(nkp_grid(d), kind=dp)
97 END DO
98
99 ALLOCATE (occupied(nkp_grid(1), nkp_grid(2), nkp_grid(3)), source=.false.)
100 DO ik = 1, nkp
101 DO d = 1, 3
102 scaled = real(nkp_grid(d), kind=dp)*(xkp(d, ik) - offset(d))
103 m = nint(scaled)
104 IF (abs(scaled - real(m, kind=dp)) > map_tolerance) RETURN
105 index(d) = modulo(m, nkp_grid(d)) + 1
106 END DO
107 IF (occupied(index(1), index(2), index(3))) RETURN
108 occupied(index(1), index(2), index(3)) = .true.
109 IF (PRESENT(grid_index)) grid_index(1:3, ik) = index
110 END DO
111 regular = .true.
112
113 IF (regular .AND. PRESENT(k_offset)) k_offset = offset
114
115 END FUNCTION regular_kpoint_grid
116
117! **************************************************************************************************
118!> \brief Find backend-supported FFT dimensions that preserve the original lattice grid.
119!> Sample the padded transform at integer multiples of the original cell indices.
120!> \param n original grid dimensions
121!> \param nfft supported integer multiples, if available
122!> \param allow_arbitrary use exact dimensions when supported by the active FFT backend
123!> \return whether all dimensions have a supported multiple
124! **************************************************************************************************
125 LOGICAL FUNCTION lattice_fft_shape(n, nfft, allow_arbitrary) RESULT(compatible)
126
127 INTEGER, DIMENSION(3), INTENT(IN) :: n
128 INTEGER, DIMENSION(3), INTENT(OUT) :: nfft
129 LOGICAL, INTENT(IN), OPTIONAL :: allow_arbitrary
130
131 INTEGER :: attempt, d, radix_length
132
133 compatible = .false.
134 nfft = n
135 IF (any(n <= 0)) RETURN
136 IF (PRESENT(allow_arbitrary)) THEN
137 IF (allow_arbitrary .AND. fft_supports_arbitrary_lengths()) THEN
138 ! For capable backends, the PW grid radix list is a cost heuristic, not a capability limit.
139 compatible = .true.
140 RETURN
141 END IF
142 END IF
143 DO d = 1, 3
144 DO attempt = 0, 15
145 IF (nfft(d) >= 3) THEN
146 CALL fft_radix_operations(nfft(d), radix_length, fft_radix_next)
147 IF (radix_length == nfft(d)) EXIT
148 END IF
149 nfft(d) = nfft(d) + n(d)
150 END DO
151 IF (attempt > 15) RETURN
152 END DO
153 compatible = .true.
154
155 END FUNCTION lattice_fft_shape
156
157! **************************************************************************************************
158!> \brief Transform a batch of real matrices from image cells to every supplied k point.
159!> \param values_rs real-space matrices; the last dimension enumerates image cells
160!> \param index_to_cell integer lattice vector for every image cell
161!> \param xkp reciprocal coordinates of all requested k points
162!> \param nkp_grid dimensions of the candidate regular reciprocal grid
163!> \param values_k complex matrices at every k point
164!> \param used_fft reports whether the regular-grid FFT path was used
165!> \param deriv_direction optional Cartesian k derivative (1, 2, or 3)
166!> \param hmat direct-lattice cell matrix, required for a derivative
167!> \param selected_kpoints optional indices of k points to retain in values_k
168! **************************************************************************************************
169 SUBROUTINE cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, values_k, &
170 used_fft, deriv_direction, hmat, selected_kpoints)
171
172 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: values_rs
173 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
174 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
175 INTEGER, DIMENSION(3), INTENT(IN) :: nkp_grid
176 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_k
177 LOGICAL, INTENT(OUT), OPTIONAL :: used_fft
178 INTEGER, INTENT(IN), OPTIONAL :: deriv_direction
179 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN), &
180 OPTIONAL :: hmat
181 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: selected_kpoints
182
183 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cell_factor
184 COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
185 POINTER :: fft_in, fft_out
186 INTEGER :: d, handle, i, icell, ik, ik_out, j, &
187 nout, stat
188 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: cell_index, grid_index
189 INTEGER, DIMENSION(3) :: n, nfft
190 LOGICAL :: regular
191 REAL(kind=dp) :: arg
192 REAL(kind=dp), DIMENSION(3) :: cell_vector, k_offset
193
194 CALL timeset("cell_to_k_grid_fft", handle)
195
196 IF (SIZE(index_to_cell, 1) < 3 .OR. &
197 SIZE(index_to_cell, 2) /= SIZE(values_rs, 3)) THEN
198 cpabort("Inconsistent real-space image-cell mapping in lattice FFT")
199 END IF
200 nout = SIZE(xkp, 2)
201 IF (PRESENT(selected_kpoints)) THEN
202 nout = SIZE(selected_kpoints)
203 IF (any(selected_kpoints < 1) .OR. any(selected_kpoints > SIZE(xkp, 2))) THEN
204 cpabort("Selected lattice-FFT k-point index is out of range")
205 END IF
206 END IF
207 IF (SIZE(values_k, 1) /= SIZE(values_rs, 1) .OR. &
208 SIZE(values_k, 2) /= SIZE(values_rs, 2) .OR. &
209 SIZE(values_k, 3) /= nout) THEN
210 cpabort("Inconsistent input and output matrix batches in lattice FFT")
211 END IF
212 IF (PRESENT(deriv_direction)) THEN
213 IF (.NOT. PRESENT(hmat)) THEN
214 CALL cp_abort(__location__, "Lattice-FFT derivative requested without a cell matrix")
215 END IF
216 IF (deriv_direction < 1 .OR. deriv_direction > 3) THEN
217 cpabort("Lattice-FFT derivative direction must be 1, 2, or 3")
218 END IF
219 END IF
220
221 ALLOCATE (grid_index(3, SIZE(xkp, 2)))
222 regular = regular_kpoint_grid(xkp, nkp_grid, grid_index, k_offset)
223 n = nkp_grid
224 IF (regular) regular = lattice_fft_shape(n, nfft)
225 IF (.NOT. regular) THEN
226 CALL direct_cell_to_k(values_rs, index_to_cell, xkp, values_k, deriv_direction, hmat, &
227 selected_kpoints)
228 IF (PRESENT(used_fft)) used_fft = .false.
229 DEALLOCATE (grid_index)
230 CALL timestop(handle)
231 RETURN
232 END IF
233
234 NULLIFY (fft_in, fft_out)
235 CALL fft_alloc(fft_in, nfft)
236 CALL fft_alloc(fft_out, nfft)
237 ALLOCATE (cell_index(3, SIZE(values_rs, 3)), cell_factor(SIZE(values_rs, 3)))
238 DO icell = 1, SIZE(values_rs, 3)
239 DO d = 1, 3
240 cell_index(d, icell) = modulo(index_to_cell(d, icell), n(d))*(nfft(d)/n(d)) + 1
241 END DO
242 arg = sum(k_offset*real(index_to_cell(1:3, icell), kind=dp))
243 cell_factor(icell) = exp(gaussi*twopi*arg)
244 IF (PRESENT(deriv_direction)) THEN
245 cell_vector = matmul(hmat, real(index_to_cell(1:3, icell), kind=dp))
246 cell_factor(icell) = cell_factor(icell)*gaussi*cell_vector(deriv_direction)
247 END IF
248 END DO
249
250 stat = 0
251 DO j = 1, SIZE(values_rs, 2)
252 DO i = 1, SIZE(values_rs, 1)
253 fft_in = z_zero
254 DO icell = 1, SIZE(values_rs, 3)
255 fft_in(cell_index(1, icell), cell_index(2, icell), cell_index(3, icell)) = &
256 fft_in(cell_index(1, icell), cell_index(2, icell), cell_index(3, icell)) + &
257 cell_factor(icell)*values_rs(i, j, icell)
258 END DO
259 CALL fft3d(bwfft, nfft, fft_in, fft_out, status=stat)
260 IF (stat /= 0) EXIT
261 DO ik_out = 1, nout
262 ik = ik_out
263 IF (PRESENT(selected_kpoints)) ik = selected_kpoints(ik_out)
264 values_k(i, j, ik_out) = &
265 fft_out(grid_index(1, ik), grid_index(2, ik), grid_index(3, ik))
266 END DO
267 END DO
268 IF (stat /= 0) EXIT
269 END DO
270
271 IF (stat /= 0) THEN
272 CALL direct_cell_to_k(values_rs, index_to_cell, xkp, values_k, deriv_direction, hmat, &
273 selected_kpoints)
274 regular = .false.
275 END IF
276 IF (PRESENT(used_fft)) used_fft = regular
277
278 CALL fft_dealloc(fft_in)
279 CALL fft_dealloc(fft_out)
280 DEALLOCATE (cell_factor, cell_index, grid_index)
281 CALL timestop(handle)
282
283 END SUBROUTINE cell_to_k_grid_fft
284
285! **************************************************************************************************
286!> \brief Transform reciprocal-grid data back to image-cell values.
287!> \param values_k complex matrices on the supplied reciprocal points
288!> \param xkp reciprocal coordinates, in units of reciprocal lattice vectors
289!> \param nkp_grid dimensions of the reciprocal grid
290!> \param index_to_cell integer lattice vector for every requested image cell
291!> \param values_rs reconstructed complex matrices in the requested image cells
292!> \param used_fft reports whether the regular-grid FFT path was used
293!> \param weights optional k weights; otherwise use the existing 1/number-of-points normalization
294!> \param allow_incomplete embed a unique regular-grid subset, zero-filling missing points
295! **************************************************************************************************
296 SUBROUTINE k_grid_to_cell_fft(values_k, xkp, nkp_grid, index_to_cell, values_rs, used_fft, &
297 weights, allow_incomplete)
298
299 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: values_k
300 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
301 INTEGER, DIMENSION(3), INTENT(IN) :: nkp_grid
302 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
303 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_rs
304 LOGICAL, INTENT(OUT), OPTIONAL :: used_fft
305 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: weights
306 LOGICAL, INTENT(IN), OPTIONAL :: allow_incomplete
307
308 TYPE(k_grid_to_cell_work_type) :: work
309
310 CALL k_grid_to_cell_prepare(work, xkp, nkp_grid, index_to_cell, weights, allow_incomplete)
311 CALL k_grid_to_cell_execute(work, values_k, values_rs, used_fft)
312 CALL k_grid_to_cell_release(work)
313
314 END SUBROUTINE k_grid_to_cell_fft
315
316! **************************************************************************************************
317!> \brief Prepare one inverse lattice transform for repeated matrix batches.
318!> \param work call-scoped maps, weights and buffers; no backend plan ownership
319!> \param xkp reciprocal coordinates
320!> \param nkp_grid candidate reciprocal-grid dimensions
321!> \param index_to_cell requested lattice translations
322!> \param weights final k weights, otherwise 1/number-of-points
323!> \param allow_incomplete allow missing grid points, with zero values
324! **************************************************************************************************
325 SUBROUTINE k_grid_to_cell_prepare(work, xkp, nkp_grid, index_to_cell, weights, allow_incomplete)
326
327 TYPE(k_grid_to_cell_work_type), INTENT(INOUT) :: work
328 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
329 INTEGER, DIMENSION(3), INTENT(IN) :: nkp_grid
330 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
331 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: weights
332 LOGICAL, INTENT(IN), OPTIONAL :: allow_incomplete
333
334 INTEGER :: d, icell, ncell, nkp
335 REAL(kind=dp), DIMENSION(3) :: k_offset
336
337 CALL k_grid_to_cell_release(work)
338 nkp = SIZE(xkp, 2)
339 ncell = SIZE(index_to_cell, 2)
340 IF (SIZE(xkp, 1) < 3 .OR. SIZE(index_to_cell, 1) < 3 .OR. nkp == 0) THEN
341 CALL cp_abort(__location__, "Inconsistent reciprocal-grid mapping in inverse lattice FFT")
342 END IF
343 ALLOCATE (work%xkp(3, nkp), work%cells(3, ncell), work%weight(nkp), work%grid_index(3, nkp))
344 work%xkp(:, :) = xkp(1:3, :)
345 work%cells(:, :) = index_to_cell(1:3, :)
346 work%weight(:) = 1.0_dp/real(nkp, dp)
347 IF (PRESENT(weights)) THEN
348 cpassert(SIZE(weights) == nkp)
349 work%weight(:) = weights
350 END IF
351 work%regular = regular_kpoint_grid(xkp, nkp_grid, work%grid_index, k_offset, allow_incomplete)
352 IF (work%regular) work%regular = lattice_fft_shape(nkp_grid, work%nfft, allow_incomplete)
353 IF (.NOT. work%regular) RETURN
354
355 CALL fft_alloc(work%fft_in, work%nfft)
356 CALL fft_alloc(work%fft_out, work%nfft)
357 ALLOCATE (work%cell_index(3, ncell), work%cell_phase(ncell))
358 DO icell = 1, ncell
359 DO d = 1, 3
360 work%cell_index(d, icell) = modulo(index_to_cell(d, icell), nkp_grid(d))* &
361 (work%nfft(d)/nkp_grid(d)) + 1
362 END DO
363 work%cell_phase(icell) = exp(-gaussi*twopi*sum(k_offset*real(index_to_cell(1:3, icell), dp)))
364 END DO
365
366 END SUBROUTINE k_grid_to_cell_prepare
367
368! **************************************************************************************************
369!> \brief Execute a matrix batch using the prepared inverse lattice transform.
370!> Call outside OpenMP worker regions: fft_tools owns its shared pool and planner.
371!> \param work prepared transform; scratch is reused, never shared between concurrent callers
372!> \param values_k matrices at the prepared reciprocal points
373!> \param values_rs matrices at the prepared lattice translations
374!> \param used_fft whether this batch used FFT, including execution failure fallback
375! **************************************************************************************************
376 SUBROUTINE k_grid_to_cell_execute(work, values_k, values_rs, used_fft)
377
378 TYPE(k_grid_to_cell_work_type), INTENT(INOUT) :: work
379 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: values_k
380 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_rs
381 LOGICAL, INTENT(OUT), OPTIONAL :: used_fft
382
383 INTEGER :: handle, i, icell, ik, j, stat
384
385 CALL timeset("k_grid_to_cell_fft", handle)
386 cpassert(ALLOCATED(work%xkp))
387 IF (SIZE(values_k, 3) /= SIZE(work%xkp, 2) .OR. &
388 SIZE(values_rs, 1) /= SIZE(values_k, 1) .OR. &
389 SIZE(values_rs, 2) /= SIZE(values_k, 2) .OR. &
390 SIZE(values_rs, 3) /= SIZE(work%cells, 2)) THEN
391 CALL cp_abort(__location__, "Inconsistent matrix batches in inverse lattice FFT")
392 END IF
393
394 stat = 0
395 IF (work%regular) THEN
396 DO j = 1, SIZE(values_k, 2)
397 DO i = 1, SIZE(values_k, 1)
398 work%fft_in = z_zero
399 DO ik = 1, SIZE(values_k, 3)
400 work%fft_in(work%grid_index(1, ik), work%grid_index(2, ik), work%grid_index(3, ik)) = &
401 work%weight(ik)*values_k(i, j, ik)
402 END DO
403 CALL fft3d(fwfft, work%nfft, work%fft_in, work%fft_out, status=stat)
404 IF (stat /= 0) EXIT
405 DO icell = 1, SIZE(values_rs, 3)
406 ! fft_tools applies 1/PRODUCT(nfft); weights are already final.
407 values_rs(i, j, icell) = real(product(work%nfft), dp)*work%cell_phase(icell)* &
408 work%fft_out(work%cell_index(1, icell), work%cell_index(2, icell), &
409 work%cell_index(3, icell))
410 END DO
411 END DO
412 IF (stat /= 0) EXIT
413 END DO
414 END IF
415 IF (stat /= 0) work%regular = .false.
416 IF (.NOT. work%regular) THEN
417 CALL direct_k_to_cell(values_k, work%xkp, work%cells, values_rs, work%weight)
418 END IF
419 IF (PRESENT(used_fft)) used_fft = work%regular
420 CALL timestop(handle)
421
422 END SUBROUTINE k_grid_to_cell_execute
423
424! **************************************************************************************************
425!> \brief Release call-scoped inverse lattice maps and FFT buffers.
426!> \param work ...
427! **************************************************************************************************
428 SUBROUTINE k_grid_to_cell_release(work)
429
430 TYPE(k_grid_to_cell_work_type), INTENT(INOUT) :: work
431
432 IF (ASSOCIATED(work%fft_in)) CALL fft_dealloc(work%fft_in)
433 IF (ASSOCIATED(work%fft_out)) CALL fft_dealloc(work%fft_out)
434 IF (ALLOCATED(work%xkp)) DEALLOCATE (work%xkp, work%cells, work%weight, work%grid_index)
435 IF (ALLOCATED(work%cell_index)) DEALLOCATE (work%cell_index, work%cell_phase)
436 work%regular = .false.
437
438 END SUBROUTINE k_grid_to_cell_release
439
440! **************************************************************************************************
441!> \brief Direct inverse phase sum for arbitrary reciprocal-point lists.
442!> \param values_k ...
443!> \param xkp ...
444!> \param index_to_cell ...
445!> \param values_rs ...
446!> \param weights optional k weights, otherwise 1/number-of-points
447! **************************************************************************************************
448 SUBROUTINE direct_k_to_cell(values_k, xkp, index_to_cell, values_rs, weights)
449
450 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: values_k
451 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
452 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
453 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_rs
454 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: weights
455
456 COMPLEX(KIND=dp) :: factor
457 INTEGER :: icell, ik
458 REAL(kind=dp) :: arg, normalization
459
460 values_rs = z_zero
461 normalization = 1.0_dp/real(SIZE(xkp, 2), kind=dp)
462 DO icell = 1, SIZE(index_to_cell, 2)
463 DO ik = 1, SIZE(xkp, 2)
464 IF (PRESENT(weights)) normalization = weights(ik)
465 arg = sum(xkp(1:3, ik)*real(index_to_cell(1:3, icell), kind=dp))
466 factor = normalization*exp(-gaussi*twopi*arg)
467 values_rs(:, :, icell) = values_rs(:, :, icell) + factor*values_k(:, :, ik)
468 END DO
469 END DO
470
471 END SUBROUTINE direct_k_to_cell
472
473! **************************************************************************************************
474!> \brief Direct phase-sum fallback for arbitrary reciprocal-point lists.
475!> \param values_rs ...
476!> \param index_to_cell ...
477!> \param xkp ...
478!> \param values_k ...
479!> \param deriv_direction ...
480!> \param hmat ...
481!> \param selected_kpoints ...
482! **************************************************************************************************
483 SUBROUTINE direct_cell_to_k(values_rs, index_to_cell, xkp, values_k, &
484 deriv_direction, hmat, selected_kpoints)
485
486 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: values_rs
487 INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
488 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
489 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: values_k
490 INTEGER, INTENT(IN), OPTIONAL :: deriv_direction
491 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN), &
492 OPTIONAL :: hmat
493 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: selected_kpoints
494
495 COMPLEX(KIND=dp) :: factor
496 INTEGER :: icell, ik, ik_out
497 REAL(kind=dp) :: arg
498 REAL(kind=dp), DIMENSION(3) :: cell_vector
499
500 values_k = z_zero
501 DO ik_out = 1, SIZE(values_k, 3)
502 ik = ik_out
503 IF (PRESENT(selected_kpoints)) ik = selected_kpoints(ik_out)
504 DO icell = 1, SIZE(values_rs, 3)
505 arg = sum(xkp(1:3, ik)*real(index_to_cell(1:3, icell), kind=dp))
506 factor = exp(gaussi*twopi*arg)
507 IF (PRESENT(deriv_direction)) THEN
508 cell_vector = matmul(hmat, real(index_to_cell(1:3, icell), kind=dp))
509 factor = factor*gaussi*cell_vector(deriv_direction)
510 END IF
511 values_k(:, :, ik_out) = values_k(:, :, ik_out) + factor*values_rs(:, :, icell)
512 END DO
513 END DO
514
515 END SUBROUTINE direct_cell_to_k
516
517END MODULE kpoint_lattice_fft
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
fft_tools::fft3d
Definition fft_tools.F:162
fft_lib
Definition fft_lib.F:7
fft_lib::fft_supports_arbitrary_lengths
logical function, public fft_supports_arbitrary_lengths()
Query whether the active FFT backend supports arbitrary transform lengths.
Definition fft_lib.F:64
fft_tools
Definition fft_tools.F:29
fft_tools::fft_radix_operations
subroutine, public fft_radix_operations(radix_in, radix_out, operation)
Determine the allowed lengths of FFT's '''.
Definition fft_tools.F:233
fft_tools::bwfft
integer, parameter, public bwfft
Definition fft_tools.F:146
fft_tools::fwfft
integer, parameter, public fwfft
Definition fft_tools.F:146
fft_tools::fft_radix_next
integer, parameter, public fft_radix_next
Definition fft_tools.F:147
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_lattice_fft
Batched lattice Fourier transforms on regular k-point grids.
Definition kpoint_lattice_fft.F:16
kpoint_lattice_fft::k_grid_to_cell_prepare
subroutine, public k_grid_to_cell_prepare(work, xkp, nkp_grid, index_to_cell, weights, allow_incomplete)
Prepare one inverse lattice transform for repeated matrix batches.
Definition kpoint_lattice_fft.F:326
kpoint_lattice_fft::regular_kpoint_grid
logical function, public regular_kpoint_grid(xkp, nkp_grid, grid_index, k_offset, allow_incomplete)
Test and map a uniformly shifted reciprocal grid or an explicitly allowed subset.
Definition kpoint_lattice_fft.F:67
kpoint_lattice_fft::k_grid_to_cell_execute
subroutine, public k_grid_to_cell_execute(work, values_k, values_rs, used_fft)
Execute a matrix batch using the prepared inverse lattice transform. Call outside OpenMP worker regio...
Definition kpoint_lattice_fft.F:377
kpoint_lattice_fft::k_grid_to_cell_release
subroutine, public k_grid_to_cell_release(work)
Release call-scoped inverse lattice maps and FFT buffers.
Definition kpoint_lattice_fft.F:429
kpoint_lattice_fft::cell_to_k_grid_fft
subroutine, public cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, values_k, used_fft, deriv_direction, hmat, selected_kpoints)
Transform a batch of real matrices from image cells to every supplied k point.
Definition kpoint_lattice_fft.F:171
kpoint_lattice_fft::lattice_fft_shape
logical function, public lattice_fft_shape(n, nfft, allow_arbitrary)
Find backend-supported FFT dimensions that preserve the original lattice grid. Sample the padded tran...
Definition kpoint_lattice_fft.F:126
kpoint_lattice_fft::k_grid_to_cell_fft
subroutine, public k_grid_to_cell_fft(values_k, xkp, nkp_grid, index_to_cell, values_rs, used_fft, weights, allow_incomplete)
Transform reciprocal-grid data back to image-cell values.
Definition kpoint_lattice_fft.F:298
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:143
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:144
kpoint_lattice_fft::k_grid_to_cell_work_type
Definition kpoint_lattice_fft.F:45