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 Implements simple k-space/R-space transformations and Gamma-point MIC reconstruction.
10 : !> \note This code is less performant/more memory consuming than the methods in kpoint_methods.F
11 : !> Use only when transformations are not the computational bottleneck.
12 : !> \par History
13 : !> 11.2025 created [Stepan Marek]
14 : !> \author Stepan Marek
15 : ! **************************************************************************************************
16 : MODULE kpoint_k_r_trafo_simple
17 : USE cell_types, ONLY: cell_type,&
18 : minimum_image_weights
19 : USE cp_cfm_types, ONLY: cp_cfm_create,&
20 : cp_cfm_type
21 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
22 : dbcsr_get_readonly_block_p,&
23 : dbcsr_p_type,&
24 : dbcsr_type,&
25 : dbcsr_type_antisymmetric,&
26 : dbcsr_type_no_symmetry
27 : USE cp_fm_struct, ONLY: cp_fm_struct_equivalent
28 : USE cp_fm_types, ONLY: cp_fm_get_info,&
29 : cp_fm_type
30 : USE kinds, ONLY: dp
31 : USE kpoint_types, ONLY: kpoint_type
32 : USE mathconstants, ONLY: twopi,&
33 : z_zero
34 : USE message_passing, ONLY: mp_comm_type
35 : USE particle_types, ONLY: particle_type
36 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
37 : neighbor_list_iterate,&
38 : neighbor_list_iterator_create,&
39 : neighbor_list_iterator_p_type,&
40 : neighbor_list_iterator_release,&
41 : neighbor_list_set_p_type
42 : #include "./base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 :
46 : PRIVATE
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_k_r_trafo_simple'
49 :
50 : PUBLIC :: replicate_rs_matrices, &
51 : rs_to_kp, &
52 : fm_gamma_to_kp, &
53 : fm_rs_to_kp, &
54 : add_kp_to_all_rs, &
55 : fm_add_kp_to_all_rs, &
56 : add_kp_to_rs
57 :
58 : CONTAINS
59 : ! **************************************************************************************************
60 : !> \brief Convert dbcsr matrices representing operators in real-space image cells to arrays
61 : !> \param rs_dbcsr_in Array of dbcsr matrices
62 : !> \param kpoint_in The kpoint environment of the source matrix (providing neighbor list and cell_to_index)
63 : !> \param rs_array_out Multidimensional array - matrices are duplicated on each MPI rank
64 : !> dim 1 : spin, dim 2 : rows, dim 3 : cols, dim 4 : image index
65 : !> \param cell_to_index_out Cell to index array for the destination array
66 : !> \date 11.2025
67 : !> \author Stepan Marek
68 : ! **************************************************************************************************
69 158 : SUBROUTINE replicate_rs_matrices(rs_dbcsr_in, kpoint_in, rs_array_out, cell_to_index_out)
70 : ! dimension 1 : spin, dimension 2 : image cell index
71 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rs_dbcsr_in
72 : TYPE(kpoint_type), POINTER :: kpoint_in
73 : REAL(kind=dp), DIMENSION(:, :, :, :) :: rs_array_out
74 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_out
75 :
76 : CHARACTER(len=*), PARAMETER :: routineN = 'replicate_rs_matrices'
77 :
78 : CHARACTER :: matrix_sym
79 : INTEGER :: handle, iatom, ispin, jatom, n_spin, &
80 : src_index
81 : INTEGER, DIMENSION(3) :: cell
82 : TYPE(mp_comm_type) :: group
83 : TYPE(neighbor_list_iterator_p_type), &
84 158 : DIMENSION(:), POINTER :: iterator
85 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
86 158 : POINTER :: sab_kp_src
87 :
88 158 : CALL timeset(routineN, handle)
89 :
90 158 : IF (SIZE(rs_dbcsr_in, 2) < 1) THEN
91 0 : CALL cp_abort(__LOCATION__, "No source image cells provided!")
92 : END IF
93 : ! Start by constructing the cell_to_index_src
94 158 : sab_kp_src => kpoint_in%sab_nl
95 : ! NOTE : The first index in matrix_s_kp is not spin index, but number of derivatives.
96 : ! But, for matrix_ks_kp, this is indeed the spin index.
97 158 : n_spin = SIZE(rs_dbcsr_in, 1)
98 158 : CALL dbcsr_get_info(rs_dbcsr_in(1, 1)%matrix, group=group, matrix_type=matrix_sym)
99 434 : DO ispin = 1, n_spin
100 276 : CALL neighbor_list_iterator_create(iterator, sab_kp_src)
101 17189 : DO WHILE (neighbor_list_iterate(iterator) == 0)
102 16913 : CALL get_iterator_info(iterator, cell=cell, iatom=iatom, jatom=jatom)
103 16913 : src_index = kpoint_in%cell_to_index(cell(1), cell(2), cell(3))
104 16913 : IF (src_index == 0) THEN
105 0 : CALL cp_abort(__LOCATION__, "Image not found in the source array.")
106 : END IF
107 : ! NOTE : Expect only specific symmetry storage relevant for kpoint calculations
108 17189 : IF (matrix_sym == dbcsr_type_no_symmetry) THEN
109 : CALL write_block_no_sym(iatom, jatom, cell, rs_dbcsr_in(ispin, src_index)%matrix, &
110 10161 : rs_array_out(ispin, :, :, :), cell_to_index_out)
111 : ELSE
112 : CALL write_block_symmetric(iatom, jatom, cell, rs_dbcsr_in(ispin, src_index)%matrix, &
113 : rs_array_out(ispin, :, :, :), cell_to_index_out, &
114 12102 : MERGE(-1.0_dp, 1.0_dp, matrix_sym == dbcsr_type_antisymmetric))
115 : END IF
116 : END DO
117 434 : CALL neighbor_list_iterator_release(iterator)
118 : END DO
119 25605570 : CALL group%sum(rs_array_out(:, :, :, :))
120 158 : CALL timestop(handle)
121 158 : END SUBROUTINE replicate_rs_matrices
122 : ! **************************************************************************************************
123 : !> \brief Expand a triangular image block and its symmetric or antisymmetric mirror.
124 : !> \param iatom first atom index
125 : !> \param jatom second atom index
126 : !> \param cell Current cell (of second atom)
127 : !> \param matrix_in Input DBCSR matrix in canonical atom order
128 : !> \param array_out Multidimensional array - dim 1 : rows, dim 2 : cols, dim 3 : rs_index
129 : !> \param cell_to_index_out Mapping of cell coords to rs_indices
130 : !> \param transpose_sign -1 for antisymmetric SOC image matrices, +1 otherwise
131 : !> \date 11.2025
132 : !> \author Stepan Marek
133 : ! **************************************************************************************************
134 6752 : SUBROUTINE write_block_symmetric(iatom, jatom, cell, matrix_in, array_out, cell_to_index_out, transpose_sign)
135 : INTEGER, INTENT(IN) :: iatom, jatom
136 : INTEGER, DIMENSION(3), INTENT(IN) :: cell
137 : TYPE(dbcsr_type), INTENT(IN) :: matrix_in
138 : REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: array_out
139 : INTEGER, DIMENSION(:, :, :), INTENT(IN), POINTER :: cell_to_index_out
140 : REAL(dp), INTENT(IN) :: transpose_sign
141 :
142 : INTEGER :: col_offset, col_size, dest_index, &
143 : dest_index_t, i, i_g, j, j_g, &
144 : row_offset, row_size
145 6752 : INTEGER, DIMENSION(:), POINTER :: col_offsets, row_offsets
146 : LOGICAL :: found
147 6752 : REAL(kind=dp), DIMENSION(:, :), POINTER :: block
148 :
149 6752 : CALL dbcsr_get_info(matrix_in, row_blk_offset=row_offsets, col_blk_offset=col_offsets)
150 6752 : IF (iatom > jatom) THEN
151 : CALL dbcsr_get_readonly_block_p(matrix_in, row=jatom, col=iatom, block=block, &
152 2333 : row_size=row_size, col_size=col_size, found=found)
153 2335 : IF (.NOT. found) RETURN
154 : ! Block found, prepare for write
155 2331 : dest_index = cell_to_index_out(-cell(1), -cell(2), -cell(3))
156 2331 : IF (dest_index == 0) CPABORT("Mirror image index not present.")
157 2331 : dest_index_t = cell_to_index_out(cell(1), cell(2), cell(3))
158 2331 : IF (dest_index_t == 0) CPABORT("Image index not present.")
159 2331 : row_offset = row_offsets(jatom)
160 2331 : col_offset = col_offsets(iatom)
161 : ELSE
162 : CALL dbcsr_get_readonly_block_p(matrix_in, row=iatom, col=jatom, block=block, &
163 4419 : row_size=row_size, col_size=col_size, found=found)
164 4419 : IF (.NOT. found) RETURN
165 : ! Block found, prepare for write
166 4417 : dest_index = cell_to_index_out(cell(1), cell(2), cell(3))
167 4417 : IF (dest_index == 0) CPABORT("Image index not present.")
168 4417 : dest_index_t = cell_to_index_out(-cell(1), -cell(2), -cell(3))
169 4417 : IF (dest_index_t == 0) CPABORT("Mirror image index not present.")
170 4417 : row_offset = row_offsets(iatom)
171 4417 : col_offset = col_offsets(jatom)
172 : END IF
173 : ! The stored block is already in canonical atom order. Only the mirror changes sign.
174 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,i_g,j_g) &
175 : !$OMP SHARED(row_size, col_size, row_offset, col_offset, array_out, &
176 6748 : !$OMP dest_index, dest_index_t, block, iatom, jatom, transpose_sign)
177 : DO i = 1, row_size
178 : i_g = i + row_offset - 1
179 : DO j = 1, col_size
180 : j_g = j + col_offset - 1
181 : array_out(i_g, j_g, dest_index) = block(i, j)
182 : IF (iatom /= jatom) array_out(j_g, i_g, dest_index_t) = transpose_sign*block(i, j)
183 : END DO
184 : END DO
185 : !$OMP END PARALLEL DO
186 6752 : END SUBROUTINE write_block_symmetric
187 : ! **************************************************************************************************
188 : !> \brief Write a single block from the dbcsr matrix to correct place, assuming no symmetry dbcsr
189 : !> \param iatom first atom index
190 : !> \param jatom second atom index
191 : !> \param cell Current cell (of second atom)
192 : !> \param matrix_in DBCSR matrix input, symmetric assumed (A_(mu nu)^(R) = A_(nu mu)^(-R))
193 : !> \param array_out Multidimensional array - dim 1 : rows, dim 2 : cols, dim 3 : rs_index
194 : !> \param cell_to_index_out Mapping of cell coords to rs_indices
195 : !> \date 11.2025
196 : !> \author Stepan Marek
197 : ! **************************************************************************************************
198 10161 : SUBROUTINE write_block_no_sym(iatom, jatom, cell, matrix_in, array_out, cell_to_index_out)
199 : INTEGER, INTENT(IN) :: iatom, jatom
200 : INTEGER, DIMENSION(3), INTENT(IN) :: cell
201 : TYPE(dbcsr_type), INTENT(IN) :: matrix_in
202 : REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: array_out
203 : INTEGER, DIMENSION(:, :, :), INTENT(IN), POINTER :: cell_to_index_out
204 :
205 : INTEGER :: col_offset, col_size, dest_index, i, &
206 : i_g, j, j_g, row_offset, row_size
207 10161 : INTEGER, DIMENSION(:), POINTER :: col_offsets, row_offsets
208 : LOGICAL :: found
209 10161 : REAL(kind=dp), DIMENSION(:, :), POINTER :: block
210 :
211 10161 : dest_index = cell_to_index_out(cell(1), cell(2), cell(3))
212 0 : IF (dest_index == 0) CPABORT("Image index not present.")
213 10161 : CALL dbcsr_get_info(matrix_in, row_blk_offset=row_offsets, col_blk_offset=col_offsets)
214 10161 : row_offset = row_offsets(iatom)
215 10161 : col_offset = col_offsets(jatom)
216 : CALL dbcsr_get_readonly_block_p(matrix_in, row=iatom, col=jatom, block=block, found=found, &
217 10161 : row_size=row_size, col_size=col_size)
218 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i, j, i_g, j_g) &
219 10161 : !$OMP SHARED(row_size, col_size, row_offset, col_offset, dest_index, block, array_out)
220 : DO i = 1, row_size
221 : i_g = i + row_offset - 1
222 : DO j = 1, col_size
223 : j_g = j + col_offset - 1
224 : array_out(i_g, j_g, dest_index) = block(i, j)
225 : END DO
226 : END DO
227 : !$OMP END PARALLEL DO
228 10161 : END SUBROUTINE write_block_no_sym
229 : ! **************************************************************************************************
230 : !> \brief Integrate RS matrices (stored as Fortran array) into a kpoint matrix at given kp
231 : !> \param rs_real Multidimensional array of real parts of the matrix, dim 1 : rows, dim 2 : cols, dim 3 : image index
232 : !> \param ks_complex Target complex k-space matrix
233 : !> \param index_to_cell Gets the image cell coordinates from the rs_dbcsr index
234 : !> \param xkp Single kpoint where the transformation is evalued at
235 : !> \param deriv_direction Derivative direction - indicates along which cartesian direction to take the derivative
236 : !> \param hmat Cell size matrix, required for the derivative
237 : !> \date 11.2025
238 : !> \author Stepan Marek
239 : ! **************************************************************************************************
240 46402 : SUBROUTINE rs_to_kp(rs_real, ks_complex, index_to_cell, xkp, deriv_direction, hmat)
241 : REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: rs_real
242 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(OUT) :: ks_complex
243 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
244 : REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xkp
245 : INTEGER, INTENT(IN), OPTIONAL :: deriv_direction
246 : REAL(kind=dp), DIMENSION(3, 3), INTENT(IN), &
247 : OPTIONAL :: hmat
248 :
249 : CHARACTER(len=*), PARAMETER :: routineN = 'rs_to_kp'
250 :
251 : INTEGER :: handle, i, n_images
252 :
253 46402 : CALL timeset(routineN, handle)
254 : ! Get matrix constants
255 46402 : n_images = SIZE(rs_real, 3)
256 : ! Get the required derivatives for the deriv direction
257 46402 : IF (PRESENT(deriv_direction)) THEN
258 603 : IF (.NOT. PRESENT(hmat) .AND. deriv_direction /= 0) THEN
259 0 : CALL cp_abort(__LOCATION__, "Derivative requested but h matrix not provided!")
260 : END IF
261 : END IF
262 : ! Now, iterate over realspace and build the sum
263 2390480 : ks_complex(:, :) = CMPLX(0.0, 0.0, kind=dp)
264 503426 : DO i = 1, n_images
265 503426 : CALL add_rs_to_kp(ks_complex, rs_real(:, :, i), xkp, i, index_to_cell, deriv_direction, hmat)
266 : END DO
267 46402 : CALL timestop(handle)
268 46402 : END SUBROUTINE rs_to_kp
269 :
270 : ! **************************************************************************************************
271 : !> \brief Reconstruct a complex k-point matrix from a real Gamma matrix using one MIC phase.
272 : !> Among images within 1.0E-6 of the minimum distance, select the last candidate
273 : !> in index_to_cell order whose lattice coordinates lie between -1 and 1.
274 : !> \param cfm_kp Complex output, with the same distributed layout as fm_gamma.
275 : !> \param fm_gamma Real Gamma matrix.
276 : !> \param cell Simulation cell.
277 : !> \param particle_set Atom positions defining the minimum images.
278 : !> \param atom_from_bf Atom index of each row/column basis function.
279 : !> \param index_to_cell Candidate lattice translations, indexed (direction, image).
280 : !> \param xkp Target k-point in reciprocal-cell coordinates.
281 : ! **************************************************************************************************
282 9524 : SUBROUTINE fm_gamma_to_kp(cfm_kp, fm_gamma, cell, particle_set, atom_from_bf, index_to_cell, xkp)
283 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_kp
284 : TYPE(cp_fm_type), INTENT(IN) :: fm_gamma
285 : TYPE(cell_type), INTENT(IN), POINTER :: cell
286 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
287 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_from_bf
288 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
289 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp
290 :
291 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_gamma_to_kp'
292 :
293 : COMPLEX(KIND=dp) :: phase
294 9524 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: phases, row_phases
295 : INTEGER :: handle, i_atom, i_atom_old, i_cell, &
296 : i_mic_cell, i_row, j_atom, j_atom_old, &
297 : j_col, ncol_local, nrow_local, &
298 : num_cells
299 9524 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
300 : LOGICAL :: same_layout
301 : REAL(KIND=dp) :: arg
302 9524 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weights
303 :
304 9524 : CALL timeset(routineN, handle)
305 :
306 9524 : CPASSERT(ASSOCIATED(fm_gamma%matrix_struct))
307 9524 : CPASSERT(fm_gamma%matrix_struct%nrow_global == fm_gamma%matrix_struct%ncol_global)
308 9524 : CPASSERT(SIZE(atom_from_bf) == fm_gamma%matrix_struct%nrow_global)
309 322004 : CPASSERT(ALL(atom_from_bf >= 1) .AND. ALL(atom_from_bf <= SIZE(particle_set)))
310 9524 : CPASSERT(SIZE(index_to_cell, 1) == 3 .AND. SIZE(index_to_cell, 2) > 0)
311 9524 : IF (.NOT. ASSOCIATED(cfm_kp%local_data)) CALL cp_cfm_create(cfm_kp, fm_gamma%matrix_struct)
312 9524 : same_layout = cp_fm_struct_equivalent(cfm_kp%matrix_struct, fm_gamma%matrix_struct)
313 95240 : same_layout = same_layout .AND. ALL(cfm_kp%matrix_struct%descriptor == fm_gamma%matrix_struct%descriptor)
314 9524 : CPASSERT(same_layout)
315 : CALL cp_fm_get_info(fm_gamma, nrow_local=nrow_local, ncol_local=ncol_local, &
316 9524 : row_indices=row_indices, col_indices=col_indices)
317 :
318 9524 : num_cells = SIZE(index_to_cell, 2)
319 66668 : ALLOCATE (weights(num_cells), phases(num_cells), row_phases(nrow_local))
320 75616 : DO i_cell = 1, num_cells
321 264368 : arg = DOT_PRODUCT(REAL(index_to_cell(:, i_cell), KIND=dp), xkp)
322 75616 : phases(i_cell) = CMPLX(COS(twopi*arg), SIN(twopi*arg), KIND=dp)
323 : END DO
324 :
325 9524 : j_atom_old = 0
326 9524 : phase = CMPLX(1.0_dp, 0.0_dp, dp)
327 165764 : DO j_col = 1, ncol_local
328 156240 : j_atom = atom_from_bf(col_indices(j_col))
329 156240 : IF (j_atom /= j_atom_old) THEN
330 : ! Reuse the atom-pair phases for all local columns on the same atom.
331 : i_atom_old = 0
332 211031 : DO i_row = 1, nrow_local
333 188777 : i_atom = atom_from_bf(row_indices(i_row))
334 188777 : IF (i_atom /= i_atom_old) THEN
335 : CALL minimum_image_weights(particle_set(i_atom)%r, particle_set(j_atom)%r, &
336 41786 : cell, index_to_cell, weights)
337 : ! The weights identify the nearest images; use only the last eligible phase.
338 : i_mic_cell = FINDLOC(weights > 0.0_dp .AND. &
339 1498036 : ALL(ABS(index_to_cell) <= 1, DIM=1), .TRUE., DIM=1, BACK=.TRUE.)
340 41786 : CPASSERT(i_mic_cell > 0)
341 41786 : phase = phases(i_mic_cell)
342 : END IF
343 188777 : row_phases(i_row) = phase
344 211031 : i_atom_old = i_atom
345 : END DO
346 : j_atom_old = j_atom
347 : END IF
348 2370514 : cfm_kp%local_data(1:nrow_local, j_col) = row_phases*fm_gamma%local_data(1:nrow_local, j_col)
349 : END DO
350 :
351 9524 : CALL timestop(handle)
352 :
353 28572 : END SUBROUTINE fm_gamma_to_kp
354 :
355 : ! **************************************************************************************************
356 : !> \brief Transforms array of fm RS matrices into cfm k-space matrix, at given kpoint index
357 : !> \param cfm_kp The resulting complex fm matrix, containing the k-space matrix elements
358 : !> \param fm_rs The real space matrix array
359 : !> \param kpoints Structure containing information about number, positions, weights etc. of kpoints
360 : !> \param ikp Index of the target kpoint at which the transformation is evaluated
361 : !> \par History
362 : !> 05.2024 Created [Jan Wilhelm]
363 : !> 11.2025 Moved to general file [Stepan Marek]
364 : ! **************************************************************************************************
365 4886 : SUBROUTINE fm_rs_to_kp(cfm_kp, fm_rs, kpoints, ikp)
366 : TYPE(cp_cfm_type) :: cfm_kp
367 : TYPE(cp_fm_type), DIMENSION(:) :: fm_rs
368 : TYPE(kpoint_type), POINTER :: kpoints
369 : INTEGER :: ikp
370 :
371 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_rs_to_kp'
372 :
373 : INTEGER :: handle, img, nimages, nimages_fm_rs
374 :
375 4886 : CALL timeset(routineN, handle)
376 :
377 4886 : nimages = SIZE(kpoints%index_to_cell, 2)
378 4886 : nimages_fm_rs = SIZE(fm_rs)
379 :
380 4886 : IF (nimages /= nimages_fm_rs) CALL cp_abort(__LOCATION__, "Index to cell and provided fm "// &
381 0 : "array are inconsistent.")
382 :
383 428743 : cfm_kp%local_data(:, :) = z_zero
384 48860 : DO img = 1, nimages
385 :
386 : CALL add_rs_to_kp(cfm_kp%local_data, fm_rs(img)%local_data, kpoints%xkp(1:3, ikp), &
387 48860 : img, kpoints%index_to_cell)
388 :
389 : END DO
390 :
391 4886 : CALL timestop(handle)
392 4886 : END SUBROUTINE fm_rs_to_kp
393 : ! **************************************************************************************************
394 : !> \brief Adds a single RS array with correct phase factor to the target k-space argument
395 : !> \param ks_array_out The array representing (part) of the complex k-space matrix
396 : !> \param rs_array_in The RS input array
397 : !> \param xkp The kpoint coordinates, in units of unit cell
398 : !> \param imindex The RS index, used for index to cell determination
399 : !> \param index_to_cell Getting cell coordinates from the imindex
400 : !> \param deriv_direction The cartesian direction of the derivative
401 : !> \param hmat The unit cell dimensions
402 : !> \par History
403 : !> 05.2024 Created [Jan Wilhelm]
404 : !> 11.2025 Added k-derivative option [Shridhar Shanbhag]
405 : !> 11.2025 Moved to general file [Stepan Marek]
406 : !> \note Always produces non-symmetric matrix, given a non-symmetric RS array
407 : ! **************************************************************************************************
408 500998 : SUBROUTINE add_rs_to_kp(ks_array_out, rs_array_in, xkp, imindex, index_to_cell, deriv_direction, hmat)
409 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: ks_array_out
410 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: rs_array_in
411 : REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xkp
412 : INTEGER, INTENT(IN) :: imindex
413 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
414 : INTEGER, INTENT(IN), OPTIONAL :: deriv_direction
415 : REAL(kind=dp), DIMENSION(3, 3), INTENT(IN), &
416 : OPTIONAL :: hmat
417 :
418 : COMPLEX(kind=dp) :: deriv_factor, phase_factor
419 : INTEGER :: i, j
420 : REAL(kind=dp), DIMENSION(3) :: cell_vector
421 :
422 500998 : IF (PRESENT(deriv_direction) .AND. (.NOT. PRESENT(hmat))) THEN
423 0 : CALL cp_abort(__LOCATION__, "Deriv. direction given but no hmat provided")
424 : END IF
425 :
426 500998 : deriv_factor = CMPLX(1.0, 0.0, kind=dp)
427 500998 : IF (PRESENT(deriv_direction)) THEN
428 485840 : cell_vector = MATMUL(hmat, index_to_cell(1:3, imindex))
429 : ! Multiply by i R(direction)
430 30365 : deriv_factor = CMPLX(0.0, cell_vector(deriv_direction), kind=dp)
431 : END IF
432 :
433 : phase_factor = CMPLX(COS(twopi*SUM(xkp(:)*index_to_cell(:, imindex))), &
434 3506986 : SIN(twopi*SUM(xkp(:)*index_to_cell(:, imindex))), kind=dp)
435 :
436 4240464 : DO i = 1, SIZE(ks_array_out, 1)
437 52623646 : DO j = 1, SIZE(ks_array_out, 2)
438 : ks_array_out(i, j) = ks_array_out(i, j) + &
439 52122648 : deriv_factor*phase_factor*rs_array_in(i, j)
440 : END DO
441 : END DO
442 500998 : END SUBROUTINE add_rs_to_kp
443 : ! **************************************************************************************************
444 : !> \brief Adds given kpoint matrix to all rs matrices
445 : !> \param array_kp Input k-space matrix array
446 : !> \param array_rs Output rs arrays
447 : !> \param kpoints Structure containing information about number, positions, weights etc. of kpoints
448 : !> \param ikp Index of the target kpoint at which the transformation is evaluated
449 : !> \param index_to_cell_ext External supplied index_to_cell
450 : !> \par History
451 : !> 05.2024 Created [Jan Wilhelm]
452 : !> 11.2025 Moved to general file [Stepan Marek]
453 : !> \note Part of transform in parallelism over k-points (hence why not for a single rs cell)
454 : ! **************************************************************************************************
455 34240 : SUBROUTINE add_kp_to_all_rs(array_kp, array_rs, kpoints, ikp, index_to_cell_ext)
456 : COMPLEX(kind=dp), DIMENSION(:, :) :: array_kp
457 : REAL(kind=dp), DIMENSION(:, :, :) :: array_rs
458 : TYPE(kpoint_type), POINTER :: kpoints
459 : INTEGER :: ikp
460 : INTEGER, DIMENSION(:, :), OPTIONAL, POINTER :: index_to_cell_ext
461 :
462 : CHARACTER(len=*), PARAMETER :: routineN = 'add_kp_to_all_rs'
463 :
464 : INTEGER :: handle, img
465 : INTEGER, DIMENSION(3) :: cell
466 34240 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
467 :
468 34240 : index_to_cell => kpoints%index_to_cell
469 33280 : IF (PRESENT(index_to_cell_ext)) index_to_cell => index_to_cell_ext
470 :
471 34240 : CALL timeset(routineN, handle)
472 :
473 34240 : IF (SIZE(array_rs, 3) /= SIZE(index_to_cell, 2)) CALL cp_abort(__LOCATION__, &
474 0 : "The provided RS array and cell_to_index array are inconsistent.")
475 :
476 342400 : DO img = 1, SIZE(array_rs, 3)
477 1232640 : cell(:) = index_to_cell(:, img)
478 342400 : CALL add_kp_to_rs(array_kp, array_rs(:, :, img), cell, kpoints, ikp)
479 : END DO
480 :
481 34240 : CALL timestop(handle)
482 34240 : END SUBROUTINE add_kp_to_all_rs
483 : ! **************************************************************************************************
484 : !> \brief Adds given kpoint matrix to a single rs matrix
485 : !> \param cfm_kp The input complex fm matrix, containing the k-space matrix elements
486 : !> \param fm_rs The real space matrix array
487 : !> \param kpoints Structure containing information about number, positions, weights etc. of kpoints
488 : !> \param ikp Index of the target kpoint at which the transformation is evaluated
489 : !> \par History
490 : !> 05.2024 Created [Jan Wilhelm]
491 : !> 11.2025 Moved to general file [Stepan Marek]
492 : !> \note Part of transform in parallelism over k-points (hence why not for a single rs cell)
493 : ! **************************************************************************************************
494 6912 : SUBROUTINE fm_add_kp_to_all_rs(cfm_kp, fm_rs, kpoints, ikp)
495 : TYPE(cp_cfm_type) :: cfm_kp
496 : TYPE(cp_fm_type), DIMENSION(:) :: fm_rs
497 : TYPE(kpoint_type), POINTER :: kpoints
498 : INTEGER :: ikp
499 :
500 : CHARACTER(len=*), PARAMETER :: routineN = 'fm_add_kp_to_all_rs'
501 :
502 : INTEGER :: handle, img
503 : INTEGER, DIMENSION(3) :: cell
504 :
505 6912 : CALL timeset(routineN, handle)
506 :
507 6912 : IF (SIZE(fm_rs, 1) /= SIZE(kpoints%index_to_cell, 2)) CALL cp_abort(__LOCATION__, &
508 0 : "The provided RS array and cell_to_index array are inconsistent.")
509 :
510 69120 : DO img = 1, SIZE(fm_rs, 1)
511 248832 : cell(:) = kpoints%index_to_cell(:, img)
512 69120 : CALL add_kp_to_rs(cfm_kp%local_data, fm_rs(img)%local_data, cell, kpoints, ikp)
513 : END DO
514 6912 : CALL timestop(handle)
515 6912 : END SUBROUTINE fm_add_kp_to_all_rs
516 : ! **************************************************************************************************
517 : !> \brief Adds given kpoint matrix to a single rs matrix
518 : !> \param array_kp The resulting complex fm matrix, containing the k-space matrix elements
519 : !> \param array_rs The real space matrix array
520 : !> \param cell The image cell of the target RS matrix
521 : !> \param kpoints Structure containing information about number, positions, weights etc. of kpoints
522 : !> \param ikp Index of the target kpoint at which the transformation is evaluated
523 : !> \par History
524 : !> 05.2024 Created [Jan Wilhelm]
525 : !> 11.2025 Moved to general file [Stepan Marek]
526 : !> \note Part of transform in parallelism over k-points (hence why not for a single rs cell)
527 : ! **************************************************************************************************
528 370368 : SUBROUTINE add_kp_to_rs(array_kp, array_rs, cell, kpoints, ikp)
529 : COMPLEX(kind=dp), DIMENSION(:, :) :: array_kp
530 : REAL(kind=dp), DIMENSION(:, :) :: array_rs
531 : INTEGER, DIMENSION(3) :: cell
532 : TYPE(kpoint_type), POINTER :: kpoints
533 : INTEGER :: ikp
534 :
535 : REAL(kind=dp) :: phase
536 :
537 : ! Determine the phase
538 1481472 : phase = twopi*SUM(kpoints%xkp(:, ikp)*cell(:))
539 : ! The array is real - it is assumed that imaginary parts cancel each other, so they are not handled
540 : ! Then e^(-i phi) S(k) = Re(S(k)) cos(phi) - i sin (phi) i Im(S(k)) =
541 : ! = Re(S(k)) cos(phi) + Im(S(k)) sin(phi)
542 : array_rs(:, :) = array_rs(:, :) + (REAL(array_kp(:, :))*COS(phase) + &
543 18832176 : AIMAG(array_kp(:, :))*SIN(phase))*kpoints%wkp(ikp)
544 370368 : END SUBROUTINE add_kp_to_rs
545 : ! **************************************************************************************************
546 : !> \brief Calculates the inverse transform, assuming all kpoints are present on the given rank
547 : !> \param array_kp The input complex fm matrix, containing the k-space matrix elements
548 : !> \param array_rs The output real space matrix array
549 : !> \param cell The image cell of the target RS matrix
550 : !> \param kpoints Structure containing information about number, positions, weights etc. of kpoints
551 : !> \par History
552 : !> 05.2024 Created [Jan Wilhelm]
553 : !> 11.2025 Moved to general file [Stepan Marek]
554 : !> \note Part of transform in parallelism over k-points
555 : ! **************************************************************************************************
556 0 : SUBROUTINE kp_to_rs(array_kp, array_rs, cell, kpoints)
557 : COMPLEX(kind=dp), DIMENSION(:, :, :) :: array_kp
558 : REAL(kind=dp), DIMENSION(:, :) :: array_rs
559 : INTEGER, DIMENSION(3) :: cell
560 : TYPE(kpoint_type), POINTER :: kpoints
561 :
562 : CHARACTER(len=*), PARAMETER :: routineN = 'kp_to_rs'
563 :
564 : INTEGER :: handle, ikp
565 :
566 0 : IF (kpoints%nkp /= SIZE(array_kp, 3)) CALL cp_abort(__LOCATION__, &
567 0 : "Provided kpoints and array_kp are inconsistent.")
568 :
569 0 : CALL timeset(routineN, handle)
570 :
571 0 : array_rs(:, :) = 0.0_dp
572 0 : DO ikp = 1, kpoints%nkp
573 0 : CALL add_kp_to_rs(array_kp(:, :, ikp), array_rs, cell, kpoints, ikp)
574 : END DO
575 :
576 0 : CALL timestop(handle)
577 0 : END SUBROUTINE kp_to_rs
578 : END MODULE kpoint_k_r_trafo_simple
|
