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 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 : ! **************************************************************************************************
16 : MODULE 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 :
55 : CONTAINS
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 120 : 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 120 : LOGICAL, ALLOCATABLE, DIMENSION(:, :, :) :: occupied
80 : REAL(KIND=dp) :: offset(3), scaled
81 :
82 120 : regular = .FALSE.
83 480 : IF (SIZE(xkp, 1) < 3 .OR. ANY(nkp_grid <= 0)) RETURN
84 120 : nkp = SIZE(xkp, 2)
85 120 : require_complete = .TRUE.
86 120 : IF (PRESENT(allow_incomplete)) require_complete = .NOT. allow_incomplete
87 480 : IF (nkp < 1 .OR. nkp > PRODUCT(nkp_grid)) RETURN
88 480 : IF (require_complete .AND. PRODUCT(nkp_grid) /= nkp) RETURN
89 120 : IF (PRESENT(grid_index)) THEN
90 44 : IF (SIZE(grid_index, 1) < 3 .OR. SIZE(grid_index, 2) < nkp) RETURN
91 : END IF
92 :
93 480 : DO d = 1, 3
94 360 : scaled = MODULO(REAL(nkp_grid(d), KIND=dp)*xkp(d, 1), 1.0_dp)
95 360 : IF (ABS(scaled) < map_tolerance .OR. ABS(scaled - 1.0_dp) < map_tolerance) scaled = 0.0_dp
96 480 : offset(d) = scaled/REAL(nkp_grid(d), KIND=dp)
97 : END DO
98 :
99 600 : ALLOCATE (occupied(nkp_grid(1), nkp_grid(2), nkp_grid(3)), source=.FALSE.)
100 7966 : DO ik = 1, nkp
101 31390 : DO d = 1, 3
102 23544 : scaled = REAL(nkp_grid(d), KIND=dp)*(xkp(d, ik) - offset(d))
103 23544 : m = NINT(scaled)
104 23544 : IF (ABS(scaled - REAL(m, KIND=dp)) > map_tolerance) RETURN
105 31388 : INDEX(d) = MODULO(m, nkp_grid(d)) + 1
106 : END DO
107 7846 : IF (occupied(INDEX(1), INDEX(2), INDEX(3))) RETURN
108 7846 : occupied(INDEX(1), INDEX(2), INDEX(3)) = .TRUE.
109 10190 : IF (PRESENT(grid_index)) grid_index(1:3, ik) = index
110 : END DO
111 118 : regular = .TRUE.
112 :
113 118 : IF (regular .AND. PRESENT(k_offset)) k_offset = offset
114 :
115 120 : 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 116 : 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 116 : compatible = .FALSE.
134 116 : nfft = n
135 464 : IF (ANY(n <= 0)) RETURN
136 116 : IF (PRESENT(allow_arbitrary)) THEN
137 82 : 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 116 : compatible = .TRUE.
140 : RETURN
141 : END IF
142 : END IF
143 160 : DO d = 1, 3
144 278 : DO attempt = 0, 15
145 278 : IF (nfft(d) >= 3) THEN
146 164 : CALL fft_radix_operations(nfft(d), radix_length, FFT_RADIX_NEXT)
147 164 : IF (radix_length == nfft(d)) EXIT
148 : END IF
149 278 : nfft(d) = nfft(d) + n(d)
150 : END DO
151 40 : IF (attempt > 15) RETURN
152 : END DO
153 116 : 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 68 : SUBROUTINE cell_to_k_grid_fft(values_rs, index_to_cell, xkp, nkp_grid, values_k, &
170 34 : 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 34 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cell_factor
184 : COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
185 34 : POINTER :: fft_in, fft_out
186 : INTEGER :: d, handle, i, icell, ik, ik_out, j, &
187 : nout, stat
188 34 : 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 34 : CALL timeset("cell_to_k_grid_fft", handle)
195 :
196 34 : IF (SIZE(index_to_cell, 1) < 3 .OR. &
197 : SIZE(index_to_cell, 2) /= SIZE(values_rs, 3)) THEN
198 0 : CPABORT("Inconsistent real-space image-cell mapping in lattice FFT")
199 : END IF
200 34 : nout = SIZE(xkp, 2)
201 34 : IF (PRESENT(selected_kpoints)) THEN
202 2 : nout = SIZE(selected_kpoints)
203 14 : IF (ANY(selected_kpoints < 1) .OR. ANY(selected_kpoints > SIZE(xkp, 2))) THEN
204 0 : 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 34 : SIZE(values_k, 2) /= SIZE(values_rs, 2) .OR. &
209 : SIZE(values_k, 3) /= nout) THEN
210 0 : CPABORT("Inconsistent input and output matrix batches in lattice FFT")
211 : END IF
212 34 : IF (PRESENT(deriv_direction)) THEN
213 2 : IF (.NOT. PRESENT(hmat)) THEN
214 0 : CALL cp_abort(__LOCATION__, "Lattice-FFT derivative requested without a cell matrix")
215 : END IF
216 2 : IF (deriv_direction < 1 .OR. deriv_direction > 3) THEN
217 0 : CPABORT("Lattice-FFT derivative direction must be 1, 2, or 3")
218 : END IF
219 : END IF
220 :
221 102 : ALLOCATE (grid_index(3, SIZE(xkp, 2)))
222 34 : regular = regular_kpoint_grid(xkp, nkp_grid, grid_index, k_offset)
223 34 : n = nkp_grid
224 34 : IF (regular) regular = lattice_fft_shape(n, nfft)
225 32 : IF (.NOT. regular) THEN
226 : CALL direct_cell_to_k(values_rs, index_to_cell, xkp, values_k, deriv_direction, hmat, &
227 4 : selected_kpoints)
228 2 : IF (PRESENT(used_fft)) used_fft = .FALSE.
229 2 : DEALLOCATE (grid_index)
230 2 : CALL timestop(handle)
231 : RETURN
232 : END IF
233 :
234 32 : NULLIFY (fft_in, fft_out)
235 32 : CALL fft_alloc(fft_in, nfft)
236 32 : CALL fft_alloc(fft_out, nfft)
237 160 : ALLOCATE (cell_index(3, SIZE(values_rs, 3)), cell_factor(SIZE(values_rs, 3)))
238 424 : DO icell = 1, SIZE(values_rs, 3)
239 1568 : DO d = 1, 3
240 1568 : cell_index(d, icell) = MODULO(index_to_cell(d, icell), n(d))*(nfft(d)/n(d)) + 1
241 : END DO
242 1568 : arg = SUM(k_offset*REAL(index_to_cell(1:3, icell), KIND=dp))
243 392 : cell_factor(icell) = EXP(gaussi*twopi*arg)
244 424 : IF (PRESENT(deriv_direction)) THEN
245 992 : cell_vector = MATMUL(hmat, REAL(index_to_cell(1:3, icell), KIND=dp))
246 62 : cell_factor(icell) = cell_factor(icell)*gaussi*cell_vector(deriv_direction)
247 : END IF
248 : END DO
249 :
250 32 : stat = 0
251 184 : DO j = 1, SIZE(values_rs, 2)
252 2614 : DO i = 1, SIZE(values_rs, 1)
253 2462 : fft_in = z_zero
254 9362 : 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 9362 : cell_factor(icell)*values_rs(i, j, icell)
258 : END DO
259 2462 : CALL fft3d(BWFFT, nfft, fft_in, fft_out, status=stat)
260 2462 : IF (stat /= 0) EXIT
261 9790 : DO ik_out = 1, nout
262 7176 : ik = ik_out
263 7176 : IF (PRESENT(selected_kpoints)) ik = selected_kpoints(ik_out)
264 : values_k(i, j, ik_out) = &
265 9638 : fft_out(grid_index(1, ik), grid_index(2, ik), grid_index(3, ik))
266 : END DO
267 : END DO
268 184 : IF (stat /= 0) EXIT
269 : END DO
270 :
271 32 : IF (stat /= 0) THEN
272 : CALL direct_cell_to_k(values_rs, index_to_cell, xkp, values_k, deriv_direction, hmat, &
273 0 : selected_kpoints)
274 0 : regular = .FALSE.
275 : END IF
276 32 : IF (PRESENT(used_fft)) used_fft = regular
277 :
278 32 : CALL fft_dealloc(fft_in)
279 32 : CALL fft_dealloc(fft_out)
280 32 : DEALLOCATE (cell_factor, cell_index, grid_index)
281 32 : CALL timestop(handle)
282 :
283 68 : 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 8 : SUBROUTINE k_grid_to_cell_fft(values_k, xkp, nkp_grid, index_to_cell, values_rs, used_fft, &
297 4 : 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 16 : TYPE(k_grid_to_cell_work_type) :: work
309 :
310 6 : CALL k_grid_to_cell_prepare(work, xkp, nkp_grid, index_to_cell, weights, allow_incomplete)
311 4 : CALL k_grid_to_cell_execute(work, values_k, values_rs, used_fft)
312 4 : CALL k_grid_to_cell_release(work)
313 :
314 4 : 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 10 : 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 10 : CALL k_grid_to_cell_release(work)
338 10 : nkp = SIZE(xkp, 2)
339 10 : ncell = SIZE(index_to_cell, 2)
340 10 : IF (SIZE(xkp, 1) < 3 .OR. SIZE(index_to_cell, 1) < 3 .OR. nkp == 0) THEN
341 0 : CALL cp_abort(__LOCATION__, "Inconsistent reciprocal-grid mapping in inverse lattice FFT")
342 : END IF
343 90 : ALLOCATE (work%xkp(3, nkp), work%cells(3, ncell), work%weight(nkp), work%grid_index(3, nkp))
344 610 : work%xkp(:, :) = xkp(1:3, :)
345 1218 : work%cells(:, :) = index_to_cell(1:3, :)
346 160 : work%weight(:) = 1.0_dp/REAL(nkp, dp)
347 10 : IF (PRESENT(weights)) THEN
348 8 : CPASSERT(SIZE(weights) == nkp)
349 126 : work%weight(:) = weights
350 : END IF
351 10 : work%regular = regular_kpoint_grid(xkp, nkp_grid, work%grid_index, k_offset, allow_incomplete)
352 10 : IF (work%regular) work%regular = lattice_fft_shape(nkp_grid, work%nfft, allow_incomplete)
353 10 : IF (.NOT. work%regular) RETURN
354 :
355 10 : CALL fft_alloc(work%fft_in, work%nfft)
356 10 : CALL fft_alloc(work%fft_out, work%nfft)
357 40 : ALLOCATE (work%cell_index(3, ncell), work%cell_phase(ncell))
358 312 : DO icell = 1, ncell
359 1208 : DO d = 1, 3
360 : work%cell_index(d, icell) = MODULO(index_to_cell(d, icell), nkp_grid(d))* &
361 1208 : (work%nfft(d)/nkp_grid(d)) + 1
362 : END DO
363 1218 : 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 172 : 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 172 : CALL timeset("k_grid_to_cell_fft", handle)
386 172 : 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 172 : SIZE(values_rs, 2) /= SIZE(values_k, 2) .OR. &
390 : SIZE(values_rs, 3) /= SIZE(work%cells, 2)) THEN
391 0 : CALL cp_abort(__LOCATION__, "Inconsistent matrix batches in inverse lattice FFT")
392 : END IF
393 :
394 172 : stat = 0
395 172 : IF (work%regular) THEN
396 364 : DO j = 1, SIZE(values_k, 2)
397 5424 : DO i = 1, SIZE(values_k, 1)
398 444720 : work%fft_in = z_zero
399 171636 : 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 171636 : work%weight(ik)*values_k(i, j, ik)
402 : END DO
403 5232 : CALL fft3d(FWFFT, work%nfft, work%fft_in, work%fft_out, status=stat)
404 5232 : IF (stat /= 0) EXIT
405 224460 : 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 881376 : work%cell_index(3, icell))
410 : END DO
411 : END DO
412 364 : IF (stat /= 0) EXIT
413 : END DO
414 : END IF
415 172 : IF (stat /= 0) work%regular = .FALSE.
416 172 : IF (.NOT. work%regular) THEN
417 0 : CALL direct_k_to_cell(values_k, work%xkp, work%cells, values_rs, work%weight)
418 : END IF
419 172 : IF (PRESENT(used_fft)) used_fft = work%regular
420 172 : CALL timestop(handle)
421 :
422 172 : 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 20 : SUBROUTINE k_grid_to_cell_release(work)
429 :
430 : TYPE(k_grid_to_cell_work_type), INTENT(INOUT) :: work
431 :
432 20 : IF (ASSOCIATED(work%fft_in)) CALL fft_dealloc(work%fft_in)
433 20 : IF (ASSOCIATED(work%fft_out)) CALL fft_dealloc(work%fft_out)
434 20 : IF (ALLOCATED(work%xkp)) DEALLOCATE (work%xkp, work%cells, work%weight, work%grid_index)
435 20 : IF (ALLOCATED(work%cell_index)) DEALLOCATE (work%cell_index, work%cell_phase)
436 20 : work%regular = .FALSE.
437 :
438 20 : 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 0 : 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 0 : values_rs = z_zero
461 0 : normalization = 1.0_dp/REAL(SIZE(xkp, 2), KIND=dp)
462 0 : DO icell = 1, SIZE(index_to_cell, 2)
463 0 : DO ik = 1, SIZE(xkp, 2)
464 0 : IF (PRESENT(weights)) normalization = weights(ik)
465 0 : arg = SUM(xkp(1:3, ik)*REAL(index_to_cell(1:3, icell), KIND=dp))
466 0 : factor = normalization*EXP(-gaussi*twopi*arg)
467 0 : values_rs(:, :, icell) = values_rs(:, :, icell) + factor*values_k(:, :, ik)
468 : END DO
469 : END DO
470 :
471 0 : 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 4 : SUBROUTINE direct_cell_to_k(values_rs, index_to_cell, xkp, values_k, &
484 2 : 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 1282 : values_k = z_zero
501 130 : DO ik_out = 1, SIZE(values_k, 3)
502 128 : ik = ik_out
503 128 : IF (PRESENT(selected_kpoints)) ik = selected_kpoints(ik_out)
504 4098 : DO icell = 1, SIZE(values_rs, 3)
505 15872 : arg = SUM(xkp(1:3, ik)*REAL(index_to_cell(1:3, icell), KIND=dp))
506 3968 : factor = EXP(gaussi*twopi*arg)
507 3968 : IF (PRESENT(deriv_direction)) THEN
508 0 : cell_vector = MATMUL(hmat, REAL(index_to_cell(1:3, icell), KIND=dp))
509 0 : factor = factor*gaussi*cell_vector(deriv_direction)
510 : END IF
511 39808 : values_k(:, :, ik_out) = values_k(:, :, ik_out) + factor*values_rs(:, :, icell)
512 : END DO
513 : END DO
514 :
515 2 : END SUBROUTINE direct_cell_to_k
516 :
517 0 : END MODULE kpoint_lattice_fft
|
