Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 24 : 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 24 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
73 24 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
74 :
75 24 : CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
76 72 : ALLOCATE (nsgfs(SIZE(particle_set)))
77 24 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
78 :
79 24 : IF (PRESENT(atom_list)) THEN
80 24 : natoms = SIZE(atom_list)
81 24 : nao = 0
82 :
83 152 : DO iatom = 1, natoms
84 152 : nao = nao + nsgfs(atom_list(iatom))
85 : END DO
86 : ELSE
87 0 : nao = SUM(nsgfs)
88 : END IF
89 :
90 24 : DEALLOCATE (nsgfs)
91 24 : 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 0 : 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 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
118 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
119 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
120 :
121 0 : CALL timeset(routineN, handle)
122 :
123 0 : CPASSERT(ASSOCIATED(matrix))
124 0 : CPASSERT(ASSOCIATED(subsys))
125 :
126 0 : CALL cp_fm_get_info(fm, nrow_global=nrows, ncol_global=ncols)
127 :
128 0 : CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
129 :
130 0 : natoms_row = SIZE(atomlist_row)
131 0 : natoms_col = SIZE(atomlist_col)
132 0 : nparticles = SIZE(particle_set)
133 :
134 0 : ALLOCATE (nsgfs(nparticles))
135 0 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
136 :
137 0 : ALLOCATE (fm_block(nrows, ncols))
138 0 : CALL cp_fm_get_submatrix(fm, fm_block)
139 :
140 0 : first_sgf_col = 1
141 0 : DO iatom_col = 1, natoms_col
142 : first_sgf_row = 1
143 0 : 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 0 : block=sm_block, found=found)
146 0 : 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 0 : nrows = nsgfs(atomlist_row(iatom_row))
151 0 : ncols = nsgfs(atomlist_col(iatom_col))
152 0 : DO icol = 1, ncols
153 0 : DO irow = 1, nrows
154 0 : 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 0 : first_sgf_row = first_sgf_row + nsgfs(atomlist_row(iatom_row))
160 : END DO
161 0 : first_sgf_col = first_sgf_col + nsgfs(atomlist_col(iatom_col))
162 : END DO
163 :
164 0 : DEALLOCATE (fm_block)
165 0 : DEALLOCATE (nsgfs)
166 :
167 0 : CALL timestop(handle)
168 0 : 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 96 : 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 96 : REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
208 : TYPE(mp_para_env_type), POINTER :: para_env
209 96 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
210 96 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
211 :
212 96 : CALL timeset(routineN, handle)
213 :
214 96 : CPASSERT(ASSOCIATED(matrix))
215 96 : CPASSERT(ASSOCIATED(subsys))
216 :
217 96 : CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
218 :
219 96 : natoms_row = SIZE(atomlist_row)
220 96 : natoms_col = SIZE(atomlist_col)
221 96 : nparticles = SIZE(particle_set)
222 :
223 288 : ALLOCATE (nsgfs(nparticles))
224 96 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
225 :
226 96 : 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 768 : CPASSERT(SUM(nsgfs(atomlist_row(:))) == nrows_fm)
230 640 : CPASSERT(SUM(nsgfs(atomlist_col(:))) == ncols_fm)
231 : END IF
232 :
233 384 : ALLOCATE (r2d(nrows_fm, ncols_fm))
234 5248 : r2d(:, :) = 0.0_dp
235 :
236 : offset_sgf_col = 0
237 640 : DO iatom_col = 1, natoms_col
238 : offset_sgf_row = 0
239 :
240 5152 : DO iatom_row = 1, natoms_row
241 4608 : IF (atomlist_row(iatom_row) <= atomlist_col(iatom_col)) THEN
242 2520 : IF (do_upper_diag) THEN
243 : CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_row(iatom_row), col=atomlist_col(iatom_col), &
244 2400 : block=sm_block, found=found)
245 : END IF
246 : ELSE
247 2088 : IF (do_lower) THEN
248 : CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_col(iatom_col), col=atomlist_row(iatom_row), &
249 2016 : block=sm_block, found=found)
250 : END IF
251 : END IF
252 :
253 4608 : IF (found) THEN
254 1594 : IF (atomlist_row(iatom_row) <= atomlist_col(iatom_col)) THEN
255 905 : IF (do_upper_diag) THEN
256 1690 : DO icol = nsgfs(atomlist_col(iatom_col)), 1, -1
257 2535 : DO irow = nsgfs(atomlist_row(iatom_row)), 1, -1
258 1690 : 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 689 : IF (do_lower) THEN
264 1306 : DO icol = nsgfs(atomlist_col(iatom_col)), 1, -1
265 1959 : DO irow = nsgfs(atomlist_row(iatom_row)), 1, -1
266 1306 : 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 5152 : offset_sgf_row = offset_sgf_row + nsgfs(atomlist_row(iatom_row))
274 : END DO
275 640 : offset_sgf_col = offset_sgf_col + nsgfs(atomlist_col(iatom_col))
276 : END DO
277 :
278 96 : CALL mpi_comm_global%sum(r2d)
279 :
280 96 : CALL cp_fm_set_submatrix(fm, r2d)
281 :
282 96 : DEALLOCATE (r2d)
283 96 : DEALLOCATE (nsgfs)
284 :
285 96 : CALL timestop(handle)
286 192 : 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 : 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 204 : IF (ABS(rep) <= 2) &
346 216 : CALL dbcsr_add(matrix_cells_raw(rep)%matrix, mat_nosym(ic)%matrix, 1.0_dp, 1.0_dp)
347 : END DO
348 :
349 12 : IF (direction_axis >= 0) THEN
350 :
351 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(1)%matrix, fm_cell1, atom_list0, atom_list1, &
352 6 : subsys, mpi_comm_global, do_upper_diag=.TRUE., do_lower=.FALSE.)
353 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(-1)%matrix, fm_cell0, atom_list0, atom_list1, &
354 6 : subsys, mpi_comm_global, do_upper_diag=.FALSE., do_lower=.TRUE.)
355 :
356 : ELSE
357 :
358 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(1)%matrix, fm_cell1, atom_list0, atom_list1, &
359 6 : subsys, mpi_comm_global, do_upper_diag=.FALSE., do_lower=.TRUE.)
360 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(-1)%matrix, fm_cell0, atom_list0, atom_list1, &
361 6 : subsys, mpi_comm_global, do_upper_diag=.TRUE., do_lower=.FALSE.)
362 :
363 : END IF
364 12 : CALL cp_fm_scale_and_add(1.0_dp, fm_cell1, 1.0_dp, fm_cell0)
365 :
366 : ! symmetric matrix fm_cell0
367 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cells_raw(0)%matrix, fm_cell0, atom_list0, atom_list0, &
368 12 : subsys, mpi_comm_global, do_upper_diag=.TRUE., do_lower=.TRUE.)
369 :
370 : ! clean up
371 12 : DEALLOCATE (c_to_i, i_to_c)
372 216 : DO ic = 1, ncell
373 216 : CALL dbcsr_release(mat_nosym(ic)%matrix)
374 : END DO
375 12 : CALL dbcsr_deallocate_matrix_set(mat_nosym)
376 72 : DO rep = -2, 2
377 72 : CALL dbcsr_deallocate_matrix(matrix_cells_raw(rep)%matrix)
378 : END DO
379 12 : DEALLOCATE (matrix_cells_raw)
380 :
381 12 : CALL timestop(handle)
382 24 : END SUBROUTINE negf_copy_contact_matrix
383 :
384 : ! **************************************************************************************************
385 : !> \brief Extract part of the DBCSR matrix based on selected atoms and copy it into another DBCSR
386 : !> matrix.
387 : !> \param matrix_contact extracted DBCSR matrix
388 : !> \param matrix_device original DBCSR matrix
389 : !> \param atom_list list of selected atoms
390 : !> \param atom_map atomic map between device and contact force environments
391 : !> \param para_env parallel environment
392 : ! **************************************************************************************************
393 0 : SUBROUTINE negf_reference_contact_matrix(matrix_contact, matrix_device, atom_list, atom_map, para_env)
394 : TYPE(dbcsr_type), POINTER :: matrix_contact, matrix_device
395 : INTEGER, DIMENSION(:), INTENT(in) :: atom_list
396 : TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
397 : TYPE(mp_para_env_type), POINTER :: para_env
398 :
399 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_reference_contact_matrix'
400 :
401 : INTEGER :: handle, i1, i2, iatom_col, iatom_row, &
402 : icol, iproc, irow, max_atom, &
403 : mepos_plus1, n1, n2, natoms, offset
404 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_nelems, send_nelems
405 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: rank_contact, rank_device
406 : LOGICAL :: found, transp
407 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: rblock
408 0 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_handlers, send_handlers
409 : TYPE(negf_allocatable_rvector), ALLOCATABLE, &
410 0 : DIMENSION(:) :: recv_packed_blocks, send_packed_blocks
411 :
412 0 : CALL timeset(routineN, handle)
413 0 : mepos_plus1 = para_env%mepos + 1
414 :
415 0 : natoms = SIZE(atom_list)
416 0 : max_atom = 0
417 0 : DO iatom_row = 1, natoms
418 0 : IF (atom_map(iatom_row)%iatom > max_atom) max_atom = atom_map(iatom_row)%iatom
419 : END DO
420 :
421 : ! find out which block goes to which node
422 0 : ALLOCATE (rank_contact(max_atom, max_atom))
423 0 : ALLOCATE (rank_device(max_atom, max_atom))
424 :
425 0 : rank_contact(:, :) = 0
426 0 : rank_device(:, :) = 0
427 :
428 0 : DO iatom_col = 1, natoms
429 0 : DO iatom_row = 1, iatom_col
430 0 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
431 0 : icol = atom_map(iatom_col)%iatom
432 0 : irow = atom_map(iatom_row)%iatom
433 : ELSE
434 0 : icol = atom_map(iatom_row)%iatom
435 0 : irow = atom_map(iatom_col)%iatom
436 : END IF
437 :
438 : CALL dbcsr_get_block_p(matrix=matrix_device, &
439 : row=atom_list(iatom_row), col=atom_list(iatom_col), &
440 0 : block=rblock, found=found)
441 0 : IF (found) rank_device(irow, icol) = mepos_plus1
442 :
443 0 : CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
444 0 : IF (found) rank_contact(irow, icol) = mepos_plus1
445 : END DO
446 : END DO
447 :
448 0 : CALL para_env%sum(rank_device)
449 0 : CALL para_env%sum(rank_contact)
450 :
451 : ! compute number of packed matrix elements to send to / receive from each processor
452 0 : ALLOCATE (recv_nelems(para_env%num_pe))
453 0 : ALLOCATE (send_nelems(para_env%num_pe))
454 0 : recv_nelems(:) = 0
455 0 : send_nelems(:) = 0
456 :
457 0 : DO iatom_col = 1, natoms
458 0 : DO iatom_row = 1, iatom_col
459 0 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
460 0 : icol = atom_map(iatom_col)%iatom
461 0 : irow = atom_map(iatom_row)%iatom
462 : ELSE
463 0 : icol = atom_map(iatom_row)%iatom
464 0 : irow = atom_map(iatom_col)%iatom
465 : END IF
466 :
467 : CALL dbcsr_get_block_p(matrix=matrix_device, &
468 : row=atom_list(iatom_row), col=atom_list(iatom_col), &
469 0 : block=rblock, found=found)
470 0 : IF (found) THEN
471 0 : iproc = rank_contact(irow, icol)
472 0 : IF (iproc > 0) &
473 0 : send_nelems(iproc) = send_nelems(iproc) + SIZE(rblock)
474 : END IF
475 :
476 0 : CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
477 0 : IF (found) THEN
478 0 : iproc = rank_device(irow, icol)
479 0 : IF (iproc > 0) &
480 0 : recv_nelems(iproc) = recv_nelems(iproc) + SIZE(rblock)
481 : END IF
482 : END DO
483 : END DO
484 :
485 : ! pack blocks
486 0 : ALLOCATE (recv_packed_blocks(para_env%num_pe))
487 0 : DO iproc = 1, para_env%num_pe
488 0 : IF (iproc /= mepos_plus1 .AND. recv_nelems(iproc) > 0) &
489 0 : ALLOCATE (recv_packed_blocks(iproc)%vector(recv_nelems(iproc)))
490 : END DO
491 :
492 0 : ALLOCATE (send_packed_blocks(para_env%num_pe))
493 0 : DO iproc = 1, para_env%num_pe
494 0 : IF (send_nelems(iproc) > 0) &
495 0 : ALLOCATE (send_packed_blocks(iproc)%vector(send_nelems(iproc)))
496 : END DO
497 :
498 0 : send_nelems(:) = 0
499 0 : DO iatom_col = 1, natoms
500 0 : DO iatom_row = 1, iatom_col
501 0 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
502 0 : icol = atom_map(iatom_col)%iatom
503 0 : irow = atom_map(iatom_row)%iatom
504 0 : transp = .FALSE.
505 : ELSE
506 0 : icol = atom_map(iatom_row)%iatom
507 0 : irow = atom_map(iatom_col)%iatom
508 0 : transp = .TRUE.
509 : END IF
510 :
511 0 : iproc = rank_contact(irow, icol)
512 0 : IF (iproc > 0) THEN
513 : CALL dbcsr_get_block_p(matrix=matrix_device, &
514 : row=atom_list(iatom_row), col=atom_list(iatom_col), &
515 0 : block=rblock, found=found)
516 0 : IF (found) THEN
517 0 : offset = send_nelems(iproc)
518 0 : n1 = SIZE(rblock, 1)
519 0 : n2 = SIZE(rblock, 2)
520 :
521 0 : IF (transp) THEN
522 0 : DO i1 = 1, n1
523 0 : DO i2 = 1, n2
524 0 : send_packed_blocks(iproc)%vector(offset + i2) = rblock(i1, i2)
525 : END DO
526 0 : offset = offset + n2
527 : END DO
528 : ELSE
529 0 : DO i2 = 1, n2
530 0 : DO i1 = 1, n1
531 0 : send_packed_blocks(iproc)%vector(offset + i1) = rblock(i1, i2)
532 : END DO
533 0 : offset = offset + n1
534 : END DO
535 : END IF
536 :
537 0 : send_nelems(iproc) = offset
538 : END IF
539 : END IF
540 : END DO
541 : END DO
542 :
543 : ! send blocks
544 0 : ALLOCATE (recv_handlers(para_env%num_pe), send_handlers(para_env%num_pe))
545 :
546 0 : DO iproc = 1, para_env%num_pe
547 0 : IF (iproc /= mepos_plus1 .AND. send_nelems(iproc) > 0) THEN
548 0 : CALL para_env%isend(send_packed_blocks(iproc)%vector, iproc - 1, send_handlers(iproc), 1)
549 : END IF
550 : END DO
551 :
552 : ! receive blocks
553 0 : DO iproc = 1, para_env%num_pe
554 0 : IF (iproc /= mepos_plus1) THEN
555 0 : IF (recv_nelems(iproc) > 0) THEN
556 0 : CALL para_env%irecv(recv_packed_blocks(iproc)%vector, iproc - 1, recv_handlers(iproc), 1)
557 : END IF
558 : ELSE
559 0 : IF (ALLOCATED(send_packed_blocks(iproc)%vector)) &
560 0 : CALL MOVE_ALLOC(send_packed_blocks(iproc)%vector, recv_packed_blocks(iproc)%vector)
561 : END IF
562 : END DO
563 :
564 : ! unpack blocks
565 0 : DO iproc = 1, para_env%num_pe
566 0 : IF (iproc /= mepos_plus1 .AND. recv_nelems(iproc) > 0) &
567 0 : CALL recv_handlers(iproc)%wait()
568 : END DO
569 :
570 0 : recv_nelems(:) = 0
571 0 : DO iatom_col = 1, natoms
572 0 : DO iatom_row = 1, iatom_col
573 0 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
574 0 : icol = atom_map(iatom_col)%iatom
575 0 : irow = atom_map(iatom_row)%iatom
576 : ELSE
577 0 : icol = atom_map(iatom_row)%iatom
578 0 : irow = atom_map(iatom_col)%iatom
579 : END IF
580 :
581 0 : iproc = rank_device(irow, icol)
582 0 : IF (iproc > 0) THEN
583 0 : CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
584 :
585 0 : IF (found) THEN
586 0 : offset = recv_nelems(iproc)
587 0 : n1 = SIZE(rblock, 1)
588 0 : n2 = SIZE(rblock, 2)
589 :
590 0 : DO i2 = 1, n2
591 0 : DO i1 = 1, n1
592 0 : rblock(i1, i2) = recv_packed_blocks(iproc)%vector(offset + i1)
593 : END DO
594 0 : offset = offset + n1
595 : END DO
596 :
597 0 : recv_nelems(iproc) = offset
598 : END IF
599 : END IF
600 : END DO
601 : END DO
602 :
603 0 : DO iproc = 1, para_env%num_pe
604 0 : IF (iproc /= mepos_plus1 .AND. send_nelems(iproc) > 0) &
605 0 : CALL send_handlers(iproc)%wait()
606 : END DO
607 :
608 : ! release memory
609 0 : DEALLOCATE (recv_handlers, send_handlers)
610 :
611 0 : DO iproc = para_env%num_pe, 1, -1
612 0 : IF (ALLOCATED(send_packed_blocks(iproc)%vector)) &
613 0 : DEALLOCATE (send_packed_blocks(iproc)%vector)
614 : END DO
615 0 : DEALLOCATE (send_packed_blocks)
616 :
617 0 : DO iproc = para_env%num_pe, 1, -1
618 0 : IF (ALLOCATED(recv_packed_blocks(iproc)%vector)) &
619 0 : DEALLOCATE (recv_packed_blocks(iproc)%vector)
620 : END DO
621 0 : DEALLOCATE (recv_packed_blocks)
622 :
623 0 : DEALLOCATE (rank_contact, rank_device)
624 0 : CALL timestop(handle)
625 0 : END SUBROUTINE negf_reference_contact_matrix
626 :
627 : ! **************************************************************************************************
628 : !> \brief Invert cell_to_index mapping between unit cells and DBCSR matrix images.
629 : !> \param cell_to_index mapping: unit_cell -> image_index
630 : !> \param nimages number of images
631 : !> \param index_to_cell inverted mapping: image_index -> unit_cell
632 : !> \par History
633 : !> * 10.2017 created [Sergey Chulkov]
634 : ! **************************************************************************************************
635 8 : SUBROUTINE invert_cell_to_index(cell_to_index, nimages, index_to_cell)
636 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
637 : INTEGER, INTENT(in) :: nimages
638 : INTEGER, DIMENSION(3, nimages), INTENT(out) :: index_to_cell
639 :
640 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invert_cell_to_index'
641 :
642 : INTEGER :: handle, i1, i2, i3, image
643 : INTEGER, DIMENSION(3) :: lbounds, ubounds
644 :
645 8 : CALL timeset(routineN, handle)
646 :
647 280 : index_to_cell(:, :) = 0
648 32 : lbounds = LBOUND(cell_to_index)
649 32 : ubounds = UBOUND(cell_to_index)
650 :
651 32 : DO i3 = lbounds(3), ubounds(3) ! z
652 96 : DO i2 = lbounds(2), ubounds(2) ! y
653 272 : DO i1 = lbounds(1), ubounds(1) ! x
654 184 : image = cell_to_index(i1, i2, i3)
655 248 : IF (image > 0 .AND. image <= nimages) THEN
656 68 : index_to_cell(1, image) = i1
657 68 : index_to_cell(2, image) = i2
658 68 : index_to_cell(3, image) = i3
659 : END IF
660 : END DO
661 : END DO
662 : END DO
663 :
664 8 : CALL timestop(handle)
665 8 : END SUBROUTINE invert_cell_to_index
666 :
667 : ! **************************************************************************************************
668 : !> \brief Helper routine to obtain index of a DBCSR matrix image by its unit cell replica.
669 : !> Can be used with any usin cell.
670 : !> \param cell indices of the unit cell
671 : !> \param cell_to_index mapping: unit_cell -> image_index
672 : !> \return DBCSR matrix images
673 : !> (0 means there are no non-zero matrix elements in the image)
674 : !> \par History
675 : !> * 10.2017 created [Sergey Chulkov]
676 : ! **************************************************************************************************
677 14888 : PURE FUNCTION get_index_by_cell(cell, cell_to_index) RESULT(image)
678 : INTEGER, DIMENSION(3), INTENT(in) :: cell
679 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
680 : INTEGER :: image
681 :
682 : IF (LBOUND(cell_to_index, 1) <= cell(1) .AND. UBOUND(cell_to_index, 1) >= cell(1) .AND. &
683 : LBOUND(cell_to_index, 2) <= cell(2) .AND. UBOUND(cell_to_index, 2) >= cell(2) .AND. &
684 104216 : LBOUND(cell_to_index, 3) <= cell(3) .AND. UBOUND(cell_to_index, 3) >= cell(3)) THEN
685 :
686 14888 : image = cell_to_index(cell(1), cell(2), cell(3))
687 : ELSE
688 : image = 0
689 : END IF
690 14888 : END FUNCTION get_index_by_cell
691 :
692 : ! **************************************************************************************************
693 : !> \brief Desymmetrizes the KS or S matrices for one of spin components
694 : !> \param mat Hamiltonian or overlap matrices
695 : !> \param mat_nosym Desymmetrized Hamiltonian or overlap matrices
696 : !> \param cell_to_index Mapping of cell indices to linear RS indices
697 : !> \param index_to_cell Mapping of linear RS indices to cell indices
698 : !> \param kpoints Kpoint environment
699 : !> \par History
700 : !> * 05.2020 created [Fabian Ducry]
701 : !> * 11.2025 Modified for one spin component. [Dmitry Ryndyk]
702 : !> \author Fabian Ducry
703 : ! **************************************************************************************************
704 12 : SUBROUTINE desymmetrize_matrix(mat, mat_nosym, cell_to_index, index_to_cell, kpoints)
705 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
706 : POINTER :: mat
707 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
708 : POINTER :: mat_nosym
709 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
710 : INTENT(OUT) :: cell_to_index
711 : INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell
712 : TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
713 :
714 : CHARACTER(len=*), PARAMETER :: routineN = 'desymmetrize_matrix'
715 :
716 : INTEGER :: handle, iatom, ic, icn, icol, irow, &
717 : jatom, ncell, nomirror, nx, ny, nz
718 : INTEGER, DIMENSION(3) :: cell
719 12 : INTEGER, DIMENSION(:, :), POINTER :: i2c
720 12 : INTEGER, DIMENSION(:, :, :), POINTER :: c2i
721 : LOGICAL :: found, lwtr
722 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
723 : TYPE(neighbor_list_iterator_p_type), &
724 12 : DIMENSION(:), POINTER :: nl_iterator
725 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
726 12 : POINTER :: sab_nl
727 :
728 12 : CALL timeset(routineN, handle)
729 :
730 12 : i2c => kpoints%index_to_cell
731 12 : c2i => kpoints%cell_to_index
732 :
733 12 : ncell = SIZE(i2c, 2)
734 :
735 12 : nx = MAX(ABS(LBOUND(c2i, 1)), ABS(UBOUND(c2i, 1)))
736 12 : ny = MAX(ABS(LBOUND(c2i, 2)), ABS(UBOUND(c2i, 3)))
737 12 : nz = MAX(ABS(LBOUND(c2i, 3)), ABS(UBOUND(c2i, 3)))
738 60 : ALLOCATE (cell_to_index(-nx:nx, -ny:ny, -nz:nz))
739 : cell_to_index(LBOUND(c2i, 1):UBOUND(c2i, 1), &
740 : LBOUND(c2i, 2):UBOUND(c2i, 2), &
741 792 : LBOUND(c2i, 3):UBOUND(c2i, 3)) = c2i
742 :
743 : ! identify cells with no mirror img
744 : nomirror = 0
745 204 : DO ic = 1, ncell
746 768 : cell = i2c(:, ic)
747 192 : IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) &
748 24 : nomirror = nomirror + 1
749 : END DO
750 :
751 : ! create the mirror imgs
752 36 : ALLOCATE (index_to_cell(3, ncell + nomirror))
753 780 : index_to_cell(:, 1:ncell) = i2c
754 :
755 : nomirror = 0 ! count the imgs without mirror
756 204 : DO ic = 1, ncell
757 768 : cell = index_to_cell(:, ic)
758 204 : IF (cell_to_index(-cell(1), -cell(2), -cell(3)) == 0) THEN
759 12 : nomirror = nomirror + 1
760 48 : index_to_cell(:, ncell + nomirror) = -cell
761 12 : cell_to_index(-cell(1), -cell(2), -cell(3)) = ncell + nomirror
762 : END IF
763 : END DO
764 12 : ncell = ncell + nomirror
765 :
766 12 : CALL get_kpoint_info(kpoints, sab_nl=sab_nl)
767 : ! allocate the nonsymmetric matrices
768 12 : NULLIFY (mat_nosym)
769 12 : CALL dbcsr_allocate_matrix_set(mat_nosym, ncell)
770 216 : DO ic = 1, ncell
771 204 : ALLOCATE (mat_nosym(ic)%matrix)
772 : CALL dbcsr_create(matrix=mat_nosym(ic)%matrix, &
773 : template=mat(1)%matrix, &
774 204 : matrix_type=dbcsr_type_no_symmetry)
775 : CALL cp_dbcsr_alloc_block_from_nbl(mat_nosym(ic)%matrix, &
776 204 : sab_nl, desymmetrize=.TRUE.)
777 216 : CALL dbcsr_set(mat_nosym(ic)%matrix, 0.0_dp)
778 : END DO
779 :
780 : ! desymmetrize the matrix for real space printing
781 12 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
782 564 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
783 552 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, cell=cell)
784 :
785 552 : ic = cell_to_index(cell(1), cell(2), cell(3))
786 552 : icn = cell_to_index(-cell(1), -cell(2), -cell(3))
787 552 : CPASSERT(icn > 0)
788 :
789 552 : irow = iatom
790 552 : icol = jatom
791 552 : lwtr = .FALSE.
792 : ! always copy from the top
793 552 : IF (iatom > jatom) THEN
794 264 : irow = jatom
795 264 : icol = iatom
796 264 : lwtr = .TRUE.
797 : END IF
798 :
799 : CALL dbcsr_get_block_p(matrix=mat(ic)%matrix, &
800 552 : row=irow, col=icol, block=block, found=found)
801 552 : CPASSERT(found)
802 :
803 : ! copy to M(R) at (iatom,jatom)
804 : ! copy to M(-R) at (jatom,iatom)
805 564 : IF (lwtr) THEN
806 : CALL dbcsr_put_block(matrix=mat_nosym(ic)%matrix, &
807 792 : row=iatom, col=jatom, block=TRANSPOSE(block))
808 : CALL dbcsr_put_block(matrix=mat_nosym(icn)%matrix, &
809 264 : row=jatom, col=iatom, block=block)
810 : ELSE
811 : CALL dbcsr_put_block(matrix=mat_nosym(ic)%matrix, &
812 288 : row=iatom, col=jatom, block=block)
813 : CALL dbcsr_put_block(matrix=mat_nosym(icn)%matrix, &
814 864 : row=jatom, col=iatom, block=TRANSPOSE(block))
815 : END IF
816 : END DO
817 12 : CALL neighbor_list_iterator_release(nl_iterator)
818 :
819 216 : DO ic = 1, ncell
820 216 : CALL dbcsr_finalize(mat_nosym(ic)%matrix)
821 : END DO
822 :
823 12 : CALL timestop(handle)
824 :
825 12 : END SUBROUTINE desymmetrize_matrix
826 :
827 : END MODULE negf_matrix_utils
|
