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 Helper routines to manipulate with matrices.
10 : ! **************************************************************************************************
11 : MODULE negf_matrix_utils
12 : USE cp_dbcsr_api, ONLY: &
13 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_finalize, &
14 : dbcsr_get_block_p, dbcsr_init_p, dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_set, &
15 : dbcsr_type, dbcsr_type_no_symmetry
16 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
17 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
18 : dbcsr_deallocate_matrix_set
19 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
20 : USE cp_fm_types, ONLY: cp_fm_get_info,&
21 : cp_fm_get_submatrix,&
22 : cp_fm_set_submatrix,&
23 : cp_fm_type
24 : USE kinds, ONLY: dp
25 : USE kpoint_types, ONLY: get_kpoint_info,&
26 : kpoint_type
27 : USE message_passing, ONLY: mp_comm_type,&
28 : mp_para_env_type,&
29 : mp_request_type
30 : USE negf_alloc_types, ONLY: negf_allocatable_rvector
31 : USE negf_atom_map, ONLY: negf_atom_map_type
32 : USE particle_methods, ONLY: get_particle_set
33 : USE particle_types, ONLY: particle_type
34 : USE qs_kind_types, ONLY: qs_kind_type
35 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
36 : neighbor_list_iterate,&
37 : neighbor_list_iterator_create,&
38 : neighbor_list_iterator_p_type,&
39 : neighbor_list_iterator_release,&
40 : neighbor_list_set_p_type
41 : USE qs_subsys_types, ONLY: qs_subsys_get,&
42 : qs_subsys_type
43 : #include "./base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 : PRIVATE
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_matrix_utils'
49 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
50 :
51 : PUBLIC :: number_of_atomic_orbitals, negf_copy_fm_submat_to_dbcsr, negf_copy_sym_dbcsr_to_fm_submat
52 : PUBLIC :: negf_copy_contact_matrix, negf_reference_contact_matrix
53 : PUBLIC :: invert_cell_to_index, get_index_by_cell
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief Compute the number of atomic orbitals of the given set of atoms.
59 : !> \param subsys QuickStep subsystem
60 : !> \param atom_list list of selected atom; when absent all the atoms are taken into account
61 : !> \return number of atomic orbitals
62 : !> \par History
63 : !> * 02.2017 created [Sergey Chulkov]
64 : ! **************************************************************************************************
65 38 : FUNCTION number_of_atomic_orbitals(subsys, atom_list) RESULT(nao)
66 : TYPE(qs_subsys_type), POINTER :: subsys
67 : INTEGER, DIMENSION(:), INTENT(in), OPTIONAL :: atom_list
68 : INTEGER :: nao
69 :
70 : INTEGER :: iatom, natoms
71 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nsgfs
72 38 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
73 38 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
74 :
75 38 : CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
76 114 : ALLOCATE (nsgfs(SIZE(particle_set)))
77 38 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
78 :
79 38 : IF (PRESENT(atom_list)) THEN
80 38 : natoms = SIZE(atom_list)
81 38 : nao = 0
82 :
83 238 : DO iatom = 1, natoms
84 238 : nao = nao + nsgfs(atom_list(iatom))
85 : END DO
86 : ELSE
87 0 : nao = SUM(nsgfs)
88 : END IF
89 :
90 38 : DEALLOCATE (nsgfs)
91 38 : END FUNCTION number_of_atomic_orbitals
92 :
93 : ! **************************************************************************************************
94 : !> \brief Populate relevant blocks of the DBCSR matrix using data from a ScaLAPACK matrix.
95 : !> Irrelevant blocks of the DBCSR matrix are kept untouched.
96 : !> \param fm dense matrix to copy
97 : !> \param matrix DBCSR matrix (modified on exit)
98 : !> \param atomlist_row set of atomic indices along the 1st (row) dimension
99 : !> \param atomlist_col set of atomic indices along the 2nd (column) dimension
100 : !> \param subsys subsystem environment
101 : !> \par History
102 : !> * 02.2017 created [Sergey Chulkov]
103 : ! **************************************************************************************************
104 20 : SUBROUTINE negf_copy_fm_submat_to_dbcsr(fm, matrix, atomlist_row, atomlist_col, subsys)
105 : TYPE(cp_fm_type), INTENT(IN) :: fm
106 : TYPE(dbcsr_type), POINTER :: matrix
107 : INTEGER, DIMENSION(:), INTENT(in) :: atomlist_row, atomlist_col
108 : TYPE(qs_subsys_type), POINTER :: subsys
109 :
110 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_copy_fm_submat_to_dbcsr'
111 :
112 : INTEGER :: first_sgf_col, first_sgf_row, handle, iatom_col, iatom_row, icol, irow, &
113 : natoms_col, natoms_row, ncols, nparticles, nrows
114 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nsgfs
115 : LOGICAL :: found
116 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: fm_block
117 20 : REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
118 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
119 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
120 :
121 20 : CALL timeset(routineN, handle)
122 :
123 20 : CPASSERT(ASSOCIATED(matrix))
124 20 : CPASSERT(ASSOCIATED(subsys))
125 :
126 20 : CALL cp_fm_get_info(fm, nrow_global=nrows, ncol_global=ncols)
127 :
128 20 : CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
129 :
130 20 : natoms_row = SIZE(atomlist_row)
131 20 : natoms_col = SIZE(atomlist_col)
132 20 : nparticles = SIZE(particle_set)
133 :
134 60 : ALLOCATE (nsgfs(nparticles))
135 20 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
136 :
137 80 : ALLOCATE (fm_block(nrows, ncols))
138 20 : CALL cp_fm_get_submatrix(fm, fm_block)
139 :
140 20 : first_sgf_col = 1
141 260 : DO iatom_col = 1, natoms_col
142 : first_sgf_row = 1
143 3120 : DO iatom_row = 1, natoms_row
144 : CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_row(iatom_row), col=atomlist_col(iatom_col), &
145 2880 : block=sm_block, found=found)
146 2880 : IF (found) THEN
147 : ! the following LAPACK call violates the coding convention
148 : !CALL dlacpy('F', nsgfs(atomlist_row(iatom_row)), nsgfs(atomlist_col(iatom_col)), &
149 : ! fm_block(first_sgf_row, first_sgf_col), SIZE(fm_block, 1), sm_block(1, 1), SIZE(sm_block, 1))
150 570 : nrows = nsgfs(atomlist_row(iatom_row))
151 570 : ncols = nsgfs(atomlist_col(iatom_col))
152 1140 : DO icol = 1, ncols
153 1710 : DO irow = 1, nrows
154 1140 : sm_block(irow, icol) = fm_block(first_sgf_row + irow - 1, first_sgf_col + icol - 1)
155 : END DO
156 : END DO
157 : END IF
158 :
159 6000 : first_sgf_row = first_sgf_row + nsgfs(atomlist_row(iatom_row))
160 : END DO
161 260 : first_sgf_col = first_sgf_col + nsgfs(atomlist_col(iatom_col))
162 : END DO
163 :
164 20 : DEALLOCATE (fm_block)
165 20 : DEALLOCATE (nsgfs)
166 :
167 20 : CALL timestop(handle)
168 40 : END SUBROUTINE negf_copy_fm_submat_to_dbcsr
169 :
170 : ! **************************************************************************************************
171 : !> \brief Extract part of the DBCSR matrix based on selected atoms and copy it into a dense matrix.
172 : !> \param matrix DBCSR matrix
173 : !> \param fm dense matrix (created and initialised on exit)
174 : !> \param atomlist_row set of atomic indices along the 1st (row) dimension
175 : !> \param atomlist_col set of atomic indices along the 2nd (column) dimension
176 : !> \param subsys subsystem environment
177 : !> \param mpi_comm_global MPI communicator which was used to distribute blocks of the DBCSR matrix.
178 : !> If missed, assume that both DBCSR and ScaLapack matrices are distributed
179 : !> across the same set of processors
180 : !> \param do_upper_diag initialise upper-triangular part of the dense matrix as well as diagonal elements
181 : !> \param do_lower initialise lower-triangular part of the dense matrix
182 : !> \par History
183 : !> * 02.2017 created [Sergey Chulkov]
184 : !> \note A naive implementation that copies relevant local DBCSR blocks into a 2-D matrix,
185 : !> performs collective summation, and then distributes the result. This approach seems to be
186 : !> optimal when processors are arranged into several independent MPI subgroups due to the fact
187 : !> that every subgroup automatically holds the copy of the dense matrix at the end, so
188 : !> we can avoid the final replication stage.
189 : ! **************************************************************************************************
190 156 : SUBROUTINE negf_copy_sym_dbcsr_to_fm_submat(matrix, fm, atomlist_row, atomlist_col, subsys, &
191 : mpi_comm_global, do_upper_diag, do_lower)
192 : TYPE(dbcsr_type), POINTER :: matrix
193 : TYPE(cp_fm_type), INTENT(IN) :: fm
194 : INTEGER, DIMENSION(:), INTENT(in) :: atomlist_row, atomlist_col
195 : TYPE(qs_subsys_type), POINTER :: subsys
196 :
197 : CLASS(mp_comm_type), INTENT(in) :: mpi_comm_global
198 : LOGICAL, INTENT(in) :: do_upper_diag, do_lower
199 :
200 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_copy_sym_dbcsr_to_fm_submat'
201 :
202 : INTEGER :: handle, iatom_col, iatom_row, icol, irow, natoms_col, natoms_row, ncols_fm, &
203 : nparticles, nrows_fm, offset_sgf_col, offset_sgf_row
204 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nsgfs
205 : LOGICAL :: found
206 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: r2d
207 156 : REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
208 : TYPE(mp_para_env_type), POINTER :: para_env
209 156 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
210 156 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
211 :
212 156 : CALL timeset(routineN, handle)
213 :
214 156 : CPASSERT(ASSOCIATED(matrix))
215 156 : CPASSERT(ASSOCIATED(subsys))
216 :
217 156 : CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
218 :
219 156 : natoms_row = SIZE(atomlist_row)
220 156 : natoms_col = SIZE(atomlist_col)
221 156 : nparticles = SIZE(particle_set)
222 :
223 468 : ALLOCATE (nsgfs(nparticles))
224 156 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
225 :
226 156 : CALL cp_fm_get_info(fm, nrow_global=nrows_fm, ncol_global=ncols_fm, para_env=para_env)
227 :
228 : IF (debug_this_module) THEN
229 1356 : CPASSERT(SUM(nsgfs(atomlist_row(:))) == nrows_fm)
230 1164 : CPASSERT(SUM(nsgfs(atomlist_col(:))) == ncols_fm)
231 : END IF
232 :
233 624 : ALLOCATE (r2d(nrows_fm, ncols_fm))
234 156 : r2d(:, :) = 0.0_dp
235 :
236 156 : offset_sgf_col = 0
237 1164 : DO iatom_col = 1, natoms_col
238 : offset_sgf_row = 0
239 :
240 10416 : DO iatom_row = 1, natoms_row
241 9408 : IF (atomlist_row(iatom_row) <= atomlist_col(iatom_col)) THEN
242 5112 : IF (do_upper_diag) THEN
243 : CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_row(iatom_row), col=atomlist_col(iatom_col), &
244 4992 : block=sm_block, found=found)
245 : END IF
246 : ELSE
247 4296 : IF (do_lower) THEN
248 : CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_col(iatom_col), col=atomlist_row(iatom_row), &
249 4224 : block=sm_block, found=found)
250 : END IF
251 : END IF
252 :
253 9408 : IF (found) THEN
254 3328 : IF (atomlist_row(iatom_row) <= atomlist_col(iatom_col)) THEN
255 1868 : IF (do_upper_diag) THEN
256 4009 : DO icol = nsgfs(atomlist_col(iatom_col)), 1, -1
257 6996 : DO irow = nsgfs(atomlist_row(iatom_row)), 1, -1
258 5188 : r2d(offset_sgf_row + irow, offset_sgf_col + icol) = sm_block(irow, icol)
259 : END DO
260 : END DO
261 : END IF
262 : ELSE
263 1460 : IF (do_lower) THEN
264 3169 : DO icol = nsgfs(atomlist_col(iatom_col)), 1, -1
265 5556 : DO irow = nsgfs(atomlist_row(iatom_row)), 1, -1
266 4132 : r2d(offset_sgf_row + irow, offset_sgf_col + icol) = sm_block(icol, irow)
267 : END DO
268 : END DO
269 : END IF
270 : END IF
271 : END IF
272 :
273 10416 : offset_sgf_row = offset_sgf_row + nsgfs(atomlist_row(iatom_row))
274 : END DO
275 1164 : offset_sgf_col = offset_sgf_col + nsgfs(atomlist_col(iatom_col))
276 : END DO
277 :
278 156 : CALL mpi_comm_global%sum(r2d)
279 :
280 156 : CALL cp_fm_set_submatrix(fm, r2d)
281 :
282 156 : DEALLOCATE (r2d)
283 156 : DEALLOCATE (nsgfs)
284 :
285 156 : CALL timestop(handle)
286 312 : END SUBROUTINE negf_copy_sym_dbcsr_to_fm_submat
287 :
288 : ! **************************************************************************************************
289 : !> \brief Driver routine to extract diagonal and off-diagonal blocks from a symmetric DBCSR matrix.
290 : !> \param fm_cell0 extracted diagonal matrix block
291 : !> \param fm_cell1 extracted off-diagonal matrix block
292 : !> \param direction_axis axis towards the secondary unit cell
293 : !> \param matrix_kp set of DBCSR matrices
294 : !> \param atom_list0 list of atoms which belong to the primary contact unit cell
295 : !> \param atom_list1 list of atoms which belong to the secondary contact unit cell
296 : !> \param subsys QuickStep subsystem
297 : !> \param mpi_comm_global global MPI communicator
298 : !> \param kpoints ...
299 : !> \par History
300 : !> * 10.2017 created [Sergey Chulkov]
301 : !> * 10.2025 The subroutine is essentially modified. [Dmitry Ryndyk]
302 : ! **************************************************************************************************
303 12 : SUBROUTINE negf_copy_contact_matrix(fm_cell0, fm_cell1, direction_axis, matrix_kp, &
304 12 : atom_list0, atom_list1, subsys, mpi_comm_global, kpoints)
305 : TYPE(cp_fm_type), INTENT(IN) :: fm_cell0, fm_cell1
306 : INTEGER, INTENT(in) :: direction_axis
307 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in), &
308 : POINTER :: matrix_kp
309 : INTEGER, DIMENSION(:), INTENT(in) :: atom_list0, atom_list1
310 : TYPE(qs_subsys_type), POINTER :: subsys
311 :
312 : CLASS(mp_comm_type), INTENT(in) :: mpi_comm_global
313 : TYPE(kpoint_type), POINTER :: kpoints
314 :
315 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_copy_contact_matrix'
316 :
317 : INTEGER :: direction_axis_abs, handle, rep, ncell, ic
318 12 : TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: matrix_cells_raw
319 12 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_nosym
320 12 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: i_to_c
321 12 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: c_to_i
322 :
323 12 : CALL timeset(routineN, handle)
324 :
325 12 : CPASSERT(ASSOCIATED(subsys))
326 :
327 12 : direction_axis_abs = ABS(direction_axis)
328 :
329 12 : CALL desymmetrize_matrix(matrix_kp, mat_nosym, c_to_i, i_to_c, kpoints)
330 12 : ncell = SIZE(i_to_c, 2) ! update the number of cells
331 :
332 : ! 0 -- primary unit cell;
333 : ! +- 1 -- upper- and lower-diagonal matrices for neighbor-cell matrix elements;
334 : ! +- 2 -- for control
335 84 : ALLOCATE (matrix_cells_raw(-2:2))
336 72 : DO rep = -2, 2
337 60 : NULLIFY (matrix_cells_raw(rep)%matrix)
338 60 : CALL dbcsr_init_p(matrix_cells_raw(rep)%matrix)
339 60 : CALL dbcsr_copy(matrix_cells_raw(rep)%matrix, mat_nosym(1)%matrix)
340 72 : CALL dbcsr_set(matrix_cells_raw(rep)%matrix, 0.0_dp)
341 : END DO
342 :
343 216 : DO ic = 1, ncell
344 204 : rep = i_to_c(direction_axis_abs, ic)
345 216 : IF (ABS(rep) <= 2) THEN
346 204 : CALL dbcsr_add(matrix_cells_raw(rep)%matrix, mat_nosym(ic)%matrix, 1.0_dp, 1.0_dp)
347 : END IF
348 : END DO
349 :
350 12 : IF (direction_axis >= 0) THEN
351 :
352 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(1)%matrix, fm_cell1, atom_list0, atom_list1, &
353 6 : subsys, mpi_comm_global, do_upper_diag=.TRUE., do_lower=.FALSE.)
354 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(-1)%matrix, fm_cell0, atom_list0, atom_list1, &
355 6 : subsys, mpi_comm_global, do_upper_diag=.FALSE., do_lower=.TRUE.)
356 :
357 : ELSE
358 :
359 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(1)%matrix, fm_cell1, atom_list0, atom_list1, &
360 6 : subsys, mpi_comm_global, do_upper_diag=.FALSE., do_lower=.TRUE.)
361 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(-1)%matrix, fm_cell0, atom_list0, atom_list1, &
362 6 : subsys, mpi_comm_global, do_upper_diag=.TRUE., do_lower=.FALSE.)
363 :
364 : END IF
365 12 : CALL cp_fm_scale_and_add(1.0_dp, fm_cell1, 1.0_dp, fm_cell0)
366 :
367 : ! symmetric matrix fm_cell0
368 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(0)%matrix, fm_cell0, atom_list0, atom_list0, &
369 12 : subsys, mpi_comm_global, do_upper_diag=.TRUE., do_lower=.TRUE.)
370 :
371 : ! clean up
372 12 : DEALLOCATE (c_to_i, i_to_c)
373 216 : DO ic = 1, ncell
374 216 : CALL dbcsr_release(mat_nosym(ic)%matrix)
375 : END DO
376 12 : CALL dbcsr_deallocate_matrix_set(mat_nosym)
377 72 : DO rep = -2, 2
378 72 : CALL dbcsr_deallocate_matrix(matrix_cells_raw(rep)%matrix)
379 : END DO
380 12 : DEALLOCATE (matrix_cells_raw)
381 :
382 12 : CALL timestop(handle)
383 12 : END SUBROUTINE negf_copy_contact_matrix
384 :
385 : ! **************************************************************************************************
386 : !> \brief Extract part of the DBCSR matrix based on selected atoms and copy it into another DBCSR
387 : !> matrix.
388 : !> \param matrix_contact extracted DBCSR matrix
389 : !> \param matrix_device original DBCSR matrix
390 : !> \param atom_list list of selected atoms
391 : !> \param atom_map atomic map between device and contact force environments
392 : !> \param para_env parallel environment
393 : ! **************************************************************************************************
394 0 : SUBROUTINE negf_reference_contact_matrix(matrix_contact, matrix_device, atom_list, atom_map, para_env)
395 : TYPE(dbcsr_type), POINTER :: matrix_contact, matrix_device
396 : INTEGER, DIMENSION(:), INTENT(in) :: atom_list
397 : TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
398 : TYPE(mp_para_env_type), POINTER :: para_env
399 :
400 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_reference_contact_matrix'
401 :
402 : INTEGER :: handle, i1, i2, iatom_col, iatom_row, &
403 : icol, iproc, irow, max_atom, &
404 : mepos_plus1, n1, n2, natoms, offset
405 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_nelems, send_nelems
406 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: rank_contact, rank_device
407 : LOGICAL :: found, transp
408 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: rblock
409 0 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_handlers, send_handlers
410 : TYPE(negf_allocatable_rvector), ALLOCATABLE, &
411 0 : DIMENSION(:) :: recv_packed_blocks, send_packed_blocks
412 :
413 0 : CALL timeset(routineN, handle)
414 0 : mepos_plus1 = para_env%mepos + 1
415 :
416 0 : natoms = SIZE(atom_list)
417 0 : max_atom = 0
418 0 : DO iatom_row = 1, natoms
419 0 : IF (atom_map(iatom_row)%iatom > max_atom) max_atom = atom_map(iatom_row)%iatom
420 : END DO
421 :
422 : ! find out which block goes to which node
423 0 : ALLOCATE (rank_contact(max_atom, max_atom))
424 0 : ALLOCATE (rank_device(max_atom, max_atom))
425 :
426 0 : rank_contact(:, :) = 0
427 0 : rank_device(:, :) = 0
428 :
429 0 : DO iatom_col = 1, natoms
430 0 : DO iatom_row = 1, iatom_col
431 0 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
432 0 : icol = atom_map(iatom_col)%iatom
433 0 : irow = atom_map(iatom_row)%iatom
434 : ELSE
435 0 : icol = atom_map(iatom_row)%iatom
436 0 : irow = atom_map(iatom_col)%iatom
437 : END IF
438 :
439 : CALL dbcsr_get_block_p(matrix=matrix_device, &
440 : row=atom_list(iatom_row), col=atom_list(iatom_col), &
441 0 : block=rblock, found=found)
442 0 : IF (found) rank_device(irow, icol) = mepos_plus1
443 :
444 0 : CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
445 0 : IF (found) rank_contact(irow, icol) = mepos_plus1
446 : END DO
447 : END DO
448 :
449 0 : CALL para_env%sum(rank_device)
450 0 : CALL para_env%sum(rank_contact)
451 :
452 : ! compute number of packed matrix elements to send to / receive from each processor
453 0 : ALLOCATE (recv_nelems(para_env%num_pe))
454 0 : ALLOCATE (send_nelems(para_env%num_pe))
455 0 : recv_nelems(:) = 0
456 0 : send_nelems(:) = 0
457 :
458 0 : DO iatom_col = 1, natoms
459 0 : DO iatom_row = 1, iatom_col
460 0 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
461 0 : icol = atom_map(iatom_col)%iatom
462 0 : irow = atom_map(iatom_row)%iatom
463 : ELSE
464 0 : icol = atom_map(iatom_row)%iatom
465 0 : irow = atom_map(iatom_col)%iatom
466 : END IF
467 :
468 : CALL dbcsr_get_block_p(matrix=matrix_device, &
469 : row=atom_list(iatom_row), col=atom_list(iatom_col), &
470 0 : block=rblock, found=found)
471 0 : IF (found) THEN
472 0 : iproc = rank_contact(irow, icol)
473 0 : IF (iproc > 0) THEN
474 0 : send_nelems(iproc) = send_nelems(iproc) + SIZE(rblock)
475 : END IF
476 : END IF
477 :
478 0 : CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
479 0 : IF (found) THEN
480 0 : iproc = rank_device(irow, icol)
481 0 : IF (iproc > 0) THEN
482 0 : recv_nelems(iproc) = recv_nelems(iproc) + SIZE(rblock)
483 : END IF
484 : END IF
485 : END DO
486 : END DO
487 :
488 : ! pack blocks
489 0 : ALLOCATE (recv_packed_blocks(para_env%num_pe))
490 0 : DO iproc = 1, para_env%num_pe
491 0 : IF (iproc /= mepos_plus1 .AND. recv_nelems(iproc) > 0) THEN
492 0 : ALLOCATE (recv_packed_blocks(iproc)%vector(recv_nelems(iproc)))
493 : END IF
494 : END DO
495 :
496 0 : ALLOCATE (send_packed_blocks(para_env%num_pe))
497 0 : DO iproc = 1, para_env%num_pe
498 0 : IF (send_nelems(iproc) > 0) THEN
499 0 : ALLOCATE (send_packed_blocks(iproc)%vector(send_nelems(iproc)))
500 : END IF
501 : END DO
502 :
503 0 : send_nelems(:) = 0
504 0 : DO iatom_col = 1, natoms
505 0 : DO iatom_row = 1, iatom_col
506 0 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
507 0 : icol = atom_map(iatom_col)%iatom
508 0 : irow = atom_map(iatom_row)%iatom
509 0 : transp = .FALSE.
510 : ELSE
511 0 : icol = atom_map(iatom_row)%iatom
512 0 : irow = atom_map(iatom_col)%iatom
513 0 : transp = .TRUE.
514 : END IF
515 :
516 0 : iproc = rank_contact(irow, icol)
517 0 : IF (iproc > 0) THEN
518 : CALL dbcsr_get_block_p(matrix=matrix_device, &
519 : row=atom_list(iatom_row), col=atom_list(iatom_col), &
520 0 : block=rblock, found=found)
521 0 : IF (found) THEN
522 0 : offset = send_nelems(iproc)
523 0 : n1 = SIZE(rblock, 1)
524 0 : n2 = SIZE(rblock, 2)
525 :
526 0 : IF (transp) THEN
527 0 : DO i1 = 1, n1
528 0 : DO i2 = 1, n2
529 0 : send_packed_blocks(iproc)%vector(offset + i2) = rblock(i1, i2)
530 : END DO
531 0 : offset = offset + n2
532 : END DO
533 : ELSE
534 0 : DO i2 = 1, n2
535 0 : DO i1 = 1, n1
536 0 : send_packed_blocks(iproc)%vector(offset + i1) = rblock(i1, i2)
537 : END DO
538 0 : offset = offset + n1
539 : END DO
540 : END IF
541 :
542 0 : send_nelems(iproc) = offset
543 : END IF
544 : END IF
545 : END DO
546 : END DO
547 :
548 : ! send blocks
549 0 : ALLOCATE (recv_handlers(para_env%num_pe), send_handlers(para_env%num_pe))
550 :
551 0 : DO iproc = 1, para_env%num_pe
552 0 : IF (iproc /= mepos_plus1 .AND. send_nelems(iproc) > 0) THEN
553 0 : CALL para_env%isend(send_packed_blocks(iproc)%vector, iproc - 1, send_handlers(iproc), 1)
554 : END IF
555 : END DO
556 :
557 : ! receive blocks
558 0 : DO iproc = 1, para_env%num_pe
559 0 : IF (iproc /= mepos_plus1) THEN
560 0 : IF (recv_nelems(iproc) > 0) THEN
561 0 : CALL para_env%irecv(recv_packed_blocks(iproc)%vector, iproc - 1, recv_handlers(iproc), 1)
562 : END IF
563 : ELSE
564 0 : IF (ALLOCATED(send_packed_blocks(iproc)%vector)) THEN
565 0 : CALL MOVE_ALLOC(send_packed_blocks(iproc)%vector, recv_packed_blocks(iproc)%vector)
566 : END IF
567 : END IF
568 : END DO
569 :
570 : ! unpack blocks
571 0 : DO iproc = 1, para_env%num_pe
572 0 : IF (iproc /= mepos_plus1 .AND. recv_nelems(iproc) > 0) THEN
573 0 : CALL recv_handlers(iproc)%wait()
574 : END IF
575 : END DO
576 :
577 0 : recv_nelems(:) = 0
578 0 : DO iatom_col = 1, natoms
579 0 : DO iatom_row = 1, iatom_col
580 0 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
581 0 : icol = atom_map(iatom_col)%iatom
582 0 : irow = atom_map(iatom_row)%iatom
583 : ELSE
584 0 : icol = atom_map(iatom_row)%iatom
585 0 : irow = atom_map(iatom_col)%iatom
586 : END IF
587 :
588 0 : iproc = rank_device(irow, icol)
589 0 : IF (iproc > 0) THEN
590 0 : CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
591 :
592 0 : IF (found) THEN
593 0 : offset = recv_nelems(iproc)
594 0 : n1 = SIZE(rblock, 1)
595 0 : n2 = SIZE(rblock, 2)
596 :
597 0 : DO i2 = 1, n2
598 0 : DO i1 = 1, n1
599 0 : rblock(i1, i2) = recv_packed_blocks(iproc)%vector(offset + i1)
600 : END DO
601 0 : offset = offset + n1
602 : END DO
603 :
604 0 : recv_nelems(iproc) = offset
605 : END IF
606 : END IF
607 : END DO
608 : END DO
609 :
610 0 : DO iproc = 1, para_env%num_pe
611 0 : IF (iproc /= mepos_plus1 .AND. send_nelems(iproc) > 0) THEN
612 0 : CALL send_handlers(iproc)%wait()
613 : END IF
614 : END DO
615 :
616 : ! release memory
617 0 : DEALLOCATE (recv_handlers, send_handlers)
618 :
619 0 : DO iproc = para_env%num_pe, 1, -1
620 0 : IF (ALLOCATED(send_packed_blocks(iproc)%vector)) THEN
621 0 : DEALLOCATE (send_packed_blocks(iproc)%vector)
622 : END IF
623 : END DO
624 0 : DEALLOCATE (send_packed_blocks)
625 :
626 0 : DO iproc = para_env%num_pe, 1, -1
627 0 : IF (ALLOCATED(recv_packed_blocks(iproc)%vector)) THEN
628 0 : DEALLOCATE (recv_packed_blocks(iproc)%vector)
629 : END IF
630 : END DO
631 0 : DEALLOCATE (recv_packed_blocks)
632 :
633 0 : DEALLOCATE (rank_contact, rank_device)
634 0 : CALL timestop(handle)
635 0 : END SUBROUTINE negf_reference_contact_matrix
636 :
637 : ! **************************************************************************************************
638 : !> \brief Invert cell_to_index mapping between unit cells and DBCSR matrix images.
639 : !> \param cell_to_index mapping: unit_cell -> image_index
640 : !> \param nimages number of images
641 : !> \param index_to_cell inverted mapping: image_index -> unit_cell
642 : !> \par History
643 : !> * 10.2017 created [Sergey Chulkov]
644 : ! **************************************************************************************************
645 16 : SUBROUTINE invert_cell_to_index(cell_to_index, nimages, index_to_cell)
646 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
647 : INTEGER, INTENT(in) :: nimages
648 : INTEGER, DIMENSION(3, nimages), INTENT(out) :: index_to_cell
649 :
650 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invert_cell_to_index'
651 :
652 : INTEGER :: handle, i1, i2, i3, image
653 : INTEGER, DIMENSION(3) :: lbounds, ubounds
654 :
655 16 : CALL timeset(routineN, handle)
656 :
657 560 : index_to_cell(:, :) = 0
658 64 : lbounds = LBOUND(cell_to_index)
659 64 : ubounds = UBOUND(cell_to_index)
660 :
661 64 : DO i3 = lbounds(3), ubounds(3) ! z
662 192 : DO i2 = lbounds(2), ubounds(2) ! y
663 544 : DO i1 = lbounds(1), ubounds(1) ! x
664 368 : image = cell_to_index(i1, i2, i3)
665 496 : IF (image > 0 .AND. image <= nimages) THEN
666 136 : index_to_cell(1, image) = i1
667 136 : index_to_cell(2, image) = i2
668 136 : index_to_cell(3, image) = i3
669 : END IF
670 : END DO
671 : END DO
672 : END DO
673 :
674 16 : CALL timestop(handle)
675 16 : END SUBROUTINE invert_cell_to_index
676 :
677 : ! **************************************************************************************************
678 : !> \brief Helper routine to obtain index of a DBCSR matrix image by its unit cell replica.
679 : !> Can be used with any usin cell.
680 : !> \param cell indices of the unit cell
681 : !> \param cell_to_index mapping: unit_cell -> image_index
682 : !> \return DBCSR matrix images
683 : !> (0 means there are no non-zero matrix elements in the image)
684 : !> \par History
685 : !> * 10.2017 created [Sergey Chulkov]
686 : ! **************************************************************************************************
687 14888 : PURE FUNCTION get_index_by_cell(cell, cell_to_index) RESULT(image)
688 : INTEGER, DIMENSION(3), INTENT(in) :: cell
689 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
690 : INTEGER :: image
691 :
692 : IF (LBOUND(cell_to_index, 1) <= cell(1) .AND. UBOUND(cell_to_index, 1) >= cell(1) .AND. &
693 : LBOUND(cell_to_index, 2) <= cell(2) .AND. UBOUND(cell_to_index, 2) >= cell(2) .AND. &
694 104216 : LBOUND(cell_to_index, 3) <= cell(3) .AND. UBOUND(cell_to_index, 3) >= cell(3)) THEN
695 :
696 14888 : image = cell_to_index(cell(1), cell(2), cell(3))
697 : ELSE
698 : image = 0
699 : END IF
700 14888 : END FUNCTION get_index_by_cell
701 :
702 : ! **************************************************************************************************
703 : !> \brief Desymmetrizes the KS or S matrices for one of spin components
704 : !> \param mat Hamiltonian or overlap matrices
705 : !> \param mat_nosym Desymmetrized Hamiltonian or overlap matrices
706 : !> \param cell_to_index Mapping of cell indices to linear RS indices
707 : !> \param index_to_cell Mapping of linear RS indices to cell indices
708 : !> \param kpoints Kpoint environment
709 : !> \par History
710 : !> * 05.2020 created [Fabian Ducry]
711 : !> * 11.2025 Modified for one spin component. [Dmitry Ryndyk]
712 : !> \author Fabian Ducry
713 : ! **************************************************************************************************
714 12 : SUBROUTINE desymmetrize_matrix(mat, mat_nosym, cell_to_index, index_to_cell, kpoints)
715 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
716 : POINTER :: mat
717 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
718 : POINTER :: mat_nosym
719 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
720 : INTENT(OUT) :: cell_to_index
721 : INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell
722 : TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
723 :
724 : CHARACTER(len=*), PARAMETER :: routineN = 'desymmetrize_matrix'
725 :
726 : INTEGER :: handle, iatom, ic, icn, icol, irow, &
727 : jatom, ncell, nomirror, nx, ny, nz
728 : INTEGER, DIMENSION(3) :: cell
729 12 : INTEGER, DIMENSION(:, :), POINTER :: i2c
730 12 : INTEGER, DIMENSION(:, :, :), POINTER :: c2i
731 : LOGICAL :: found, lwtr
732 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
733 : TYPE(neighbor_list_iterator_p_type), &
734 12 : DIMENSION(:), POINTER :: nl_iterator
735 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
736 12 : POINTER :: sab_nl
737 :
738 12 : CALL timeset(routineN, handle)
739 :
740 12 : i2c => kpoints%index_to_cell
741 12 : c2i => kpoints%cell_to_index
742 :
743 12 : ncell = SIZE(i2c, 2)
744 :
745 12 : nx = MAX(ABS(LBOUND(c2i, 1)), ABS(UBOUND(c2i, 1)))
746 12 : ny = MAX(ABS(LBOUND(c2i, 2)), ABS(UBOUND(c2i, 3)))
747 12 : nz = MAX(ABS(LBOUND(c2i, 3)), ABS(UBOUND(c2i, 3)))
748 60 : ALLOCATE (cell_to_index(-nx:nx, -ny:ny, -nz:nz))
749 : cell_to_index(LBOUND(c2i, 1):UBOUND(c2i, 1), &
750 : LBOUND(c2i, 2):UBOUND(c2i, 2), &
751 792 : LBOUND(c2i, 3):UBOUND(c2i, 3)) = c2i
752 :
753 : ! identify cells with no mirror img
754 : nomirror = 0
755 204 : DO ic = 1, ncell
756 768 : cell = i2c(:, ic)
757 204 : IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) THEN
758 12 : nomirror = nomirror + 1
759 : END IF
760 : END DO
761 :
762 : ! create the mirror imgs
763 36 : ALLOCATE (index_to_cell(3, ncell + nomirror))
764 780 : index_to_cell(:, 1:ncell) = i2c
765 :
766 : nomirror = 0 ! count the imgs without mirror
767 204 : DO ic = 1, ncell
768 768 : cell = index_to_cell(:, ic)
769 204 : IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) THEN
770 12 : nomirror = nomirror + 1
771 48 : index_to_cell(:, ncell + nomirror) = -cell
772 12 : cell_to_index(-cell(1), -cell(2), -cell(3)) = ncell + nomirror
773 : END IF
774 : END DO
775 12 : ncell = ncell + nomirror
776 :
777 12 : CALL get_kpoint_info(kpoints, sab_nl=sab_nl)
778 : ! allocate the nonsymmetric matrices
779 12 : NULLIFY (mat_nosym)
780 12 : CALL dbcsr_allocate_matrix_set(mat_nosym, ncell)
781 216 : DO ic = 1, ncell
782 204 : ALLOCATE (mat_nosym(ic)%matrix)
783 : CALL dbcsr_create(matrix=mat_nosym(ic)%matrix, &
784 : template=mat(1)%matrix, &
785 204 : matrix_type=dbcsr_type_no_symmetry)
786 : CALL cp_dbcsr_alloc_block_from_nbl(mat_nosym(ic)%matrix, &
787 204 : sab_nl, desymmetrize=.TRUE.)
788 216 : CALL dbcsr_set(mat_nosym(ic)%matrix, 0.0_dp)
789 : END DO
790 :
791 : ! desymmetrize the matrix for real space printing
792 12 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
793 564 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
794 552 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
795 :
796 552 : ic = cell_to_index(cell(1), cell(2), cell(3))
797 552 : icn = cell_to_index(-cell(1), -cell(2), -cell(3))
798 552 : CPASSERT(icn > 0)
799 :
800 552 : irow = iatom
801 552 : icol = jatom
802 552 : lwtr = .FALSE.
803 : ! always copy from the top
804 552 : IF (iatom > jatom) THEN
805 264 : irow = jatom
806 264 : icol = iatom
807 264 : lwtr = .TRUE.
808 : END IF
809 :
810 : CALL dbcsr_get_block_p(matrix=mat(ic)%matrix, &
811 552 : row=irow, col=icol, block=block, found=found)
812 552 : CPASSERT(found)
813 :
814 : ! copy to M(R) at (iatom,jatom)
815 : ! copy to M(-R) at (jatom,iatom)
816 564 : IF (lwtr) THEN
817 : CALL dbcsr_put_block(matrix=mat_nosym(ic)%matrix, &
818 792 : row=iatom, col=jatom, block=TRANSPOSE(block))
819 : CALL dbcsr_put_block(matrix=mat_nosym(icn)%matrix, &
820 264 : row=jatom, col=iatom, block=block)
821 : ELSE
822 : CALL dbcsr_put_block(matrix=mat_nosym(ic)%matrix, &
823 288 : row=iatom, col=jatom, block=block)
824 : CALL dbcsr_put_block(matrix=mat_nosym(icn)%matrix, &
825 864 : row=jatom, col=iatom, block=TRANSPOSE(block))
826 : END IF
827 : END DO
828 12 : CALL neighbor_list_iterator_release(nl_iterator)
829 :
830 216 : DO ic = 1, ncell
831 216 : CALL dbcsr_finalize(mat_nosym(ic)%matrix)
832 : END DO
833 :
834 12 : CALL timestop(handle)
835 :
836 12 : END SUBROUTINE desymmetrize_matrix
837 :
838 : END MODULE negf_matrix_utils
|
