(git:374b731)
Loading...
Searching...
No Matches
negf_matrix_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 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
12MODULE negf_matrix_utils
13 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
14 USE cp_fm_types, ONLY: cp_fm_get_info,&
15 cp_fm_get_submatrix,&
16 cp_fm_set_submatrix,&
17 cp_fm_type
18 USE dbcsr_api, ONLY: &
19 dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_get_block_p, dbcsr_get_info, &
20 dbcsr_init_p, dbcsr_p_type, dbcsr_set, dbcsr_type
21 USE kinds, ONLY: dp
22 USE message_passing, ONLY: mp_comm_type,&
23 mp_para_env_type,&
24 mp_request_type
25 USE negf_alloc_types, ONLY: negf_allocatable_rvector
26 USE negf_atom_map, ONLY: negf_atom_map_type
27 USE particle_methods, ONLY: get_particle_set
28 USE particle_types, ONLY: particle_type
29 USE qs_kind_types, ONLY: qs_kind_type
30 USE qs_subsys_types, ONLY: qs_subsys_get,&
31 qs_subsys_type
32#include "./base/base_uses.f90"
33
34 IMPLICIT NONE
35 PRIVATE
36
37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_matrix_utils'
38 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .true.
39
40 PUBLIC :: number_of_atomic_orbitals, negf_copy_fm_submat_to_dbcsr, negf_copy_sym_dbcsr_to_fm_submat
41 PUBLIC :: negf_copy_contact_matrix, negf_reference_contact_matrix
42 PUBLIC :: invert_cell_to_index, get_index_by_cell
43
44CONTAINS
45
46! **************************************************************************************************
47!> \brief Compute the number of atomic orbitals of the given set of atoms.
48!> \param subsys QuickStep subsystem
49!> \param atom_list list of selected atom; when absent all the atoms are taken into account
50!> \return number of atomic orbitals
51!> \par History
52!> * 02.2017 created [Sergey Chulkov]
53! **************************************************************************************************
54 FUNCTION number_of_atomic_orbitals(subsys, atom_list) RESULT(nao)
55 TYPE(qs_subsys_type), POINTER :: subsys
56 INTEGER, DIMENSION(:), INTENT(in), OPTIONAL :: atom_list
57 INTEGER :: nao
58
59 INTEGER :: iatom, natoms
60 INTEGER, ALLOCATABLE, DIMENSION(:) :: nsgfs
61 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
62 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
63
64 CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
65 ALLOCATE (nsgfs(SIZE(particle_set)))
66 CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
67
68 IF (PRESENT(atom_list)) THEN
69 natoms = SIZE(atom_list)
70 nao = 0
71
72 DO iatom = 1, natoms
73 nao = nao + nsgfs(atom_list(iatom))
74 END DO
75 ELSE
76 nao = sum(nsgfs)
77 END IF
78
79 DEALLOCATE (nsgfs)
80 END FUNCTION number_of_atomic_orbitals
81
82! **************************************************************************************************
83!> \brief Populate relevant blocks of the DBCSR matrix using data from a ScaLAPACK matrix.
84!> Irrelevant blocks of the DBCSR matrix are kept untouched.
85!> \param fm dense matrix to copy
86!> \param matrix DBCSR matrix (modified on exit)
87!> \param atomlist_row set of atomic indices along the 1st (row) dimension
88!> \param atomlist_col set of atomic indices along the 2nd (column) dimension
89!> \param subsys subsystem environment
90!> \par History
91!> * 02.2017 created [Sergey Chulkov]
92! **************************************************************************************************
93 SUBROUTINE negf_copy_fm_submat_to_dbcsr(fm, matrix, atomlist_row, atomlist_col, subsys)
94 TYPE(cp_fm_type), INTENT(IN) :: fm
95 TYPE(dbcsr_type), POINTER :: matrix
96 INTEGER, DIMENSION(:), INTENT(in) :: atomlist_row, atomlist_col
97 TYPE(qs_subsys_type), POINTER :: subsys
98
99 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_copy_fm_submat_to_dbcsr'
100
101 INTEGER :: first_sgf_col, first_sgf_row, handle, iatom_col, iatom_row, icol, irow, &
102 natoms_col, natoms_row, ncols, nparticles, nrows
103 INTEGER, ALLOCATABLE, DIMENSION(:) :: nsgfs
104 LOGICAL :: found
105 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: fm_block
106 REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
107 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
109
110 CALL timeset(routinen, handle)
111
112 cpassert(ASSOCIATED(matrix))
113 cpassert(ASSOCIATED(subsys))
114
115 CALL cp_fm_get_info(fm, nrow_global=nrows, ncol_global=ncols)
116
117 CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
118
119 natoms_row = SIZE(atomlist_row)
120 natoms_col = SIZE(atomlist_col)
121 nparticles = SIZE(particle_set)
122
123 ALLOCATE (nsgfs(nparticles))
124 CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
125
126 ALLOCATE (fm_block(nrows, ncols))
127 CALL cp_fm_get_submatrix(fm, fm_block)
128
129 first_sgf_col = 1
130 DO iatom_col = 1, natoms_col
131 first_sgf_row = 1
132 DO iatom_row = 1, natoms_row
133 CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_row(iatom_row), col=atomlist_col(iatom_col), &
134 block=sm_block, found=found)
135 IF (found) THEN
136 ! the following LAPACK call violates the coding convention
137 !CALL dlacpy('F', nsgfs(atomlist_row(iatom_row)), nsgfs(atomlist_col(iatom_col)), &
138 ! fm_block(first_sgf_row, first_sgf_col), SIZE(fm_block, 1), sm_block(1, 1), SIZE(sm_block, 1))
139 nrows = nsgfs(atomlist_row(iatom_row))
140 ncols = nsgfs(atomlist_col(iatom_col))
141 DO icol = 1, ncols
142 DO irow = 1, nrows
143 sm_block(irow, icol) = fm_block(first_sgf_row + irow - 1, first_sgf_col + icol - 1)
144 END DO
145 END DO
146 END IF
147
148 first_sgf_row = first_sgf_row + nsgfs(atomlist_row(iatom_row))
149 END DO
150 first_sgf_col = first_sgf_col + nsgfs(atomlist_col(iatom_col))
151 END DO
152
153 DEALLOCATE (fm_block)
154 DEALLOCATE (nsgfs)
155
156 CALL timestop(handle)
157 END SUBROUTINE negf_copy_fm_submat_to_dbcsr
158
159! **************************************************************************************************
160!> \brief Extract part of the DBCSR matrix based on selected atoms and copy it into a dense matrix.
161!> \param matrix DBCSR matrix
162!> \param fm dense matrix (created and initialised on exit)
163!> \param atomlist_row set of atomic indices along the 1st (row) dimension
164!> \param atomlist_col set of atomic indices along the 2nd (column) dimension
165!> \param subsys subsystem environment
166!> \param mpi_comm_global MPI communicator which was used to distribute blocks of the DBCSR matrix.
167!> If missed, assume that both DBCSR and ScaLapack matrices are distributed
168!> across the same set of processors
169!> \param do_upper_diag initialise upper-triangular part of the dense matrix as well as diagonal elements
170!> \param do_lower initialise lower-triangular part of the dense matrix
171!> \par History
172!> * 02.2017 created [Sergey Chulkov]
173!> \note A naive implementation that copies relevant local DBCSR blocks into a 2-D matrix,
174!> performs collective summation, and then distributes the result. This approach seems to be
175!> optimal when processors are arranged into several independent MPI subgroups due to the fact
176!> that every subgroup automatically holds the copy of the dense matrix at the end, so
177!> we can avoid the final replication stage.
178! **************************************************************************************************
179 SUBROUTINE negf_copy_sym_dbcsr_to_fm_submat(matrix, fm, atomlist_row, atomlist_col, subsys, &
180 mpi_comm_global, do_upper_diag, do_lower)
181 TYPE(dbcsr_type), POINTER :: matrix
182 TYPE(cp_fm_type), INTENT(IN) :: fm
183 INTEGER, DIMENSION(:), INTENT(in) :: atomlist_row, atomlist_col
184 TYPE(qs_subsys_type), POINTER :: subsys
185
186 CLASS(mp_comm_type), INTENT(in) :: mpi_comm_global
187 LOGICAL, INTENT(in) :: do_upper_diag, do_lower
188
189 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_copy_sym_dbcsr_to_fm_submat'
190
191 INTEGER :: handle, iatom_col, iatom_row, icol, irow, natoms_col, natoms_row, ncols_fm, &
192 nparticles, nrows_fm, offset_sgf_col, offset_sgf_row
193 INTEGER, ALLOCATABLE, DIMENSION(:) :: nsgfs
194 LOGICAL :: found
195 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: r2d
196 REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
197 TYPE(mp_para_env_type), POINTER :: para_env
198 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
199 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
200
201 CALL timeset(routinen, handle)
202
203 cpassert(ASSOCIATED(matrix))
204 cpassert(ASSOCIATED(subsys))
205
206 CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
207
208 natoms_row = SIZE(atomlist_row)
209 natoms_col = SIZE(atomlist_col)
210 nparticles = SIZE(particle_set)
211
212 ALLOCATE (nsgfs(nparticles))
213 CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
214
215 CALL cp_fm_get_info(fm, nrow_global=nrows_fm, ncol_global=ncols_fm, para_env=para_env)
216
217 IF (debug_this_module) THEN
218 cpassert(sum(nsgfs(atomlist_row(:))) == nrows_fm)
219 cpassert(sum(nsgfs(atomlist_col(:))) == ncols_fm)
220 END IF
221
222 ALLOCATE (r2d(nrows_fm, ncols_fm))
223 r2d(:, :) = 0.0_dp
224
225 offset_sgf_col = 0
226 DO iatom_col = 1, natoms_col
227 offset_sgf_row = 0
228
229 DO iatom_row = 1, natoms_row
230 IF (atomlist_row(iatom_row) <= atomlist_col(iatom_col)) THEN
231 IF (do_upper_diag) THEN
232 CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_row(iatom_row), col=atomlist_col(iatom_col), &
233 block=sm_block, found=found)
234 END IF
235 ELSE
236 IF (do_lower) THEN
237 CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_col(iatom_col), col=atomlist_row(iatom_row), &
238 block=sm_block, found=found)
239 END IF
240 END IF
241
242 IF (found) THEN
243 IF (atomlist_row(iatom_row) <= atomlist_col(iatom_col)) THEN
244 IF (do_upper_diag) THEN
245 DO icol = nsgfs(atomlist_col(iatom_col)), 1, -1
246 DO irow = nsgfs(atomlist_row(iatom_row)), 1, -1
247 r2d(offset_sgf_row + irow, offset_sgf_col + icol) = sm_block(irow, icol)
248 END DO
249 END DO
250 END IF
251 ELSE
252 IF (do_lower) THEN
253 DO icol = nsgfs(atomlist_col(iatom_col)), 1, -1
254 DO irow = nsgfs(atomlist_row(iatom_row)), 1, -1
255 r2d(offset_sgf_row + irow, offset_sgf_col + icol) = sm_block(icol, irow)
256 END DO
257 END DO
258 END IF
259 END IF
260 END IF
261
262 offset_sgf_row = offset_sgf_row + nsgfs(atomlist_row(iatom_row))
263 END DO
264 offset_sgf_col = offset_sgf_col + nsgfs(atomlist_col(iatom_col))
265 END DO
266
267 CALL mpi_comm_global%sum(r2d)
268
269 CALL cp_fm_set_submatrix(fm, r2d)
270
271 DEALLOCATE (r2d)
272 DEALLOCATE (nsgfs)
273
274 CALL timestop(handle)
275 END SUBROUTINE negf_copy_sym_dbcsr_to_fm_submat
276
277! **************************************************************************************************
278!> \brief Driver routine to extract diagonal and off-diagonal blocks from a symmetric DBCSR matrix.
279!> \param fm_cell0 extracted diagonal matrix block
280!> \param fm_cell1 extracted off-diagonal matrix block
281!> \param direction_axis axis towards the secondary unit cell
282!> \param matrix_kp set of DBCSR matrices
283!> \param index_to_cell inverted mapping between unit cells and DBCSR matrix images
284!> \param atom_list0 list of atoms which belong to the primary contact unit cell
285!> \param atom_list1 list of atoms which belong to the secondary contact unit cell
286!> \param subsys QuickStep subsystem
287!> \param mpi_comm_global global MPI communicator
288!> \param is_same_cell for every atomic pair indicates whether or not both atoms are assigned to
289!> the same (0) or different (-1) unit cells (initialised when the optional
290!> argument 'matrix_ref' is given)
291!> \param matrix_ref reference DBCSR matrix
292!> \par History
293!> * 10.2017 created [Sergey Chulkov]
294! **************************************************************************************************
295 SUBROUTINE negf_copy_contact_matrix(fm_cell0, fm_cell1, direction_axis, matrix_kp, index_to_cell, &
296 atom_list0, atom_list1, subsys, mpi_comm_global, is_same_cell, matrix_ref)
297 TYPE(cp_fm_type), INTENT(IN) :: fm_cell0, fm_cell1
298 INTEGER, INTENT(in) :: direction_axis
299 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_kp
300 INTEGER, DIMENSION(:, :), INTENT(in) :: index_to_cell
301 INTEGER, DIMENSION(:), INTENT(in) :: atom_list0, atom_list1
302 TYPE(qs_subsys_type), POINTER :: subsys
303
304 CLASS(mp_comm_type), INTENT(in) :: mpi_comm_global
305 INTEGER, DIMENSION(:, :), INTENT(inout) :: is_same_cell
306 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_ref
307
308 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_copy_contact_matrix'
309
310 INTEGER :: direction_axis_abs, handle, iatom_col, &
311 iatom_row, image, natoms, nimages, &
312 phase, rep
313 LOGICAL :: found
314 REAL(kind=dp) :: error_diff, error_same
315 REAL(kind=dp), DIMENSION(:, :), POINTER :: block_dest, block_src
316 TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: matrix_cells_raw
317 TYPE(dbcsr_type), POINTER :: matrix_cell_0, matrix_cell_1, &
318 matrix_cell_minus1
319
320 CALL timeset(routinen, handle)
321
322 cpassert(ASSOCIATED(subsys))
323
324 nimages = SIZE(index_to_cell, 2)
325 direction_axis_abs = abs(direction_axis)
326
327 ! 0 -- primary unit cell;
328 ! +- 1 -- upper- and lower-diagonal matrices for the secondary unit cell;
329 ! when the distance between two atoms within the unit cell becomes bigger than
330 ! the distance between the same atoms from different cell replicas, the third
331 ! unit cell replica (+- 2) is also needed.
332 ALLOCATE (matrix_cells_raw(-2:2))
333 DO rep = -2, 2
334 NULLIFY (matrix_cells_raw(rep)%matrix)
335 CALL dbcsr_init_p(matrix_cells_raw(rep)%matrix)
336 CALL dbcsr_copy(matrix_cells_raw(rep)%matrix, matrix_kp(1)%matrix)
337 CALL dbcsr_set(matrix_cells_raw(rep)%matrix, 0.0_dp)
338 END DO
339
340 NULLIFY (matrix_cell_0, matrix_cell_1, matrix_cell_minus1)
341
342 CALL dbcsr_init_p(matrix_cell_0)
343 CALL dbcsr_copy(matrix_cell_0, matrix_kp(1)%matrix)
344 CALL dbcsr_set(matrix_cell_0, 0.0_dp)
345
346 CALL dbcsr_init_p(matrix_cell_1)
347 CALL dbcsr_copy(matrix_cell_1, matrix_kp(1)%matrix)
348 CALL dbcsr_set(matrix_cell_1, 0.0_dp)
349
350 CALL dbcsr_init_p(matrix_cell_minus1)
351 CALL dbcsr_copy(matrix_cell_minus1, matrix_kp(1)%matrix)
352 CALL dbcsr_set(matrix_cell_minus1, 0.0_dp)
353
354 DO image = 1, nimages
355 rep = index_to_cell(direction_axis_abs, image)
356
357 IF (abs(rep) <= 2) &
358 CALL dbcsr_add(matrix_cells_raw(rep)%matrix, matrix_kp(image)%matrix, 1.0_dp, 1.0_dp)
359 END DO
360
361 CALL dbcsr_get_info(matrix_cell_0, nblkrows_total=natoms)
362
363 IF (PRESENT(matrix_ref)) THEN
364 ! 0 -- atoms belong to the same cell or absent (zero) matrix block;
365 ! +1 -- atoms belong to different cells
366 is_same_cell(:, :) = 0
367
368 DO iatom_col = 1, natoms
369 DO iatom_row = 1, iatom_col
370 CALL dbcsr_get_block_p(matrix=matrix_ref, &
371 row=iatom_row, col=iatom_col, &
372 block=block_src, found=found)
373
374 IF (found) THEN
375 ! it should be much safe to rely on atomic indices (iatom / jatom) obtained using a neighbour list iterator:
376 ! phase == 1 when iatom <= jatom, and phase == -1 when iatom > jatom
377 IF (mod(iatom_col - iatom_row, 2) == 0) THEN
378 phase = 1
379 ELSE
380 phase = -1
381 END IF
382
383 CALL dbcsr_get_block_p(matrix=matrix_cells_raw(0)%matrix, &
384 row=iatom_row, col=iatom_col, &
385 block=block_dest, found=found)
386 cpassert(found)
387
388 error_same = maxval(abs(block_dest(:, :) - block_src(:, :)))
389
390 CALL dbcsr_get_block_p(matrix=matrix_cells_raw(phase)%matrix, &
391 row=iatom_row, col=iatom_col, &
392 block=block_dest, found=found)
393 cpassert(found)
394 error_diff = maxval(abs(block_dest(:, :) - block_src(:, :)))
395
396 IF (error_same <= error_diff) THEN
397 is_same_cell(iatom_row, iatom_col) = 0
398 ELSE
399 is_same_cell(iatom_row, iatom_col) = 1
400 END IF
401 END IF
402 END DO
403 END DO
404 END IF
405
406 DO iatom_col = 1, natoms
407 DO iatom_row = 1, iatom_col
408 CALL dbcsr_get_block_p(matrix=matrix_cell_0, &
409 row=iatom_row, col=iatom_col, block=block_dest, found=found)
410
411 IF (found) THEN
412 ! it should be much safe to rely on a neighbour list iterator
413 IF (mod(iatom_col - iatom_row, 2) == 0) THEN
414 phase = 1
415 ELSE
416 phase = -1
417 END IF
418 rep = phase*is_same_cell(iatom_row, iatom_col)
419
420 ! primary unit cell:
421 ! matrix(i,j) <- [0]%matrix(i,j) when i and j are from the same replica
422 ! matrix(i,j) <- [phase]%matrix(i,j) when i and j are from different replicas
423 CALL dbcsr_get_block_p(matrix=matrix_cells_raw(rep)%matrix, &
424 row=iatom_row, col=iatom_col, block=block_src, found=found)
425 cpassert(found)
426 block_dest(:, :) = block_src(:, :)
427
428 ! secondary unit cell, i <= j:
429 ! matrix(i,j) <- [phase]%matrix(i,j) when i and j are from the same replica
430 ! matrix(i,j) <- [2*phase]%matrix(i,j) when i and j are from different replicas
431 CALL dbcsr_get_block_p(matrix=matrix_cell_1, &
432 row=iatom_row, col=iatom_col, block=block_dest, found=found)
433 cpassert(found)
434 CALL dbcsr_get_block_p(matrix=matrix_cells_raw(rep + phase)%matrix, &
435 row=iatom_row, col=iatom_col, block=block_src, found=found)
436 cpassert(found)
437 block_dest(:, :) = block_src(:, :)
438
439 ! secondary unit cell, i > j:
440 ! matrix(i,j) <- [-phase]%matrix(i,j) when i and j are from the same replica
441 ! matrix(i,j) <- [-2*phase]%matrix(i,j) when i and j are from different replicas
442 CALL dbcsr_get_block_p(matrix=matrix_cell_minus1, &
443 row=iatom_row, col=iatom_col, block=block_dest, found=found)
444 cpassert(found)
445 CALL dbcsr_get_block_p(matrix=matrix_cells_raw(rep - phase)%matrix, &
446 row=iatom_row, col=iatom_col, block=block_src, found=found)
447 cpassert(found)
448 block_dest(:, :) = block_src(:, :)
449 END IF
450 END DO
451 END DO
452
453 IF (direction_axis >= 0) THEN
454 ! upper-diagonal part of fm_cell1
455 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cell_1, fm_cell1, atom_list0, atom_list1, &
456 subsys, mpi_comm_global, do_upper_diag=.true., do_lower=.false.)
457 ! lower-diagonal part of fm_cell1
458 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cell_minus1, fm_cell0, atom_list0, atom_list1, &
459 subsys, mpi_comm_global, do_upper_diag=.false., do_lower=.true.)
460 ELSE
461 ! upper-diagonal part of fm_cell1
462 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cell_minus1, fm_cell1, atom_list0, atom_list1, &
463 subsys, mpi_comm_global, do_upper_diag=.true., do_lower=.false.)
464 ! lower-diagonal part of fm_cell1
465 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cell_1, fm_cell0, atom_list0, atom_list1, &
466 subsys, mpi_comm_global, do_upper_diag=.false., do_lower=.true.)
467
468 END IF
469 CALL cp_fm_scale_and_add(1.0_dp, fm_cell1, 1.0_dp, fm_cell0)
470
471 ! symmetric matrix fm_cell0
472 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix_cell_0, fm_cell0, atom_list0, atom_list0, &
473 subsys, mpi_comm_global, do_upper_diag=.true., do_lower=.true.)
474
475 CALL dbcsr_deallocate_matrix(matrix_cell_0)
476 CALL dbcsr_deallocate_matrix(matrix_cell_1)
477 CALL dbcsr_deallocate_matrix(matrix_cell_minus1)
478
479 DO rep = -2, 2
480 CALL dbcsr_deallocate_matrix(matrix_cells_raw(rep)%matrix)
481 END DO
482 DEALLOCATE (matrix_cells_raw)
483
484 CALL timestop(handle)
485 END SUBROUTINE negf_copy_contact_matrix
486
487! **************************************************************************************************
488!> \brief Extract part of the DBCSR matrix based on selected atoms and copy it into another DBCSR
489!> matrix.
490!> \param matrix_contact extracted DBCSR matrix
491!> \param matrix_device original DBCSR matrix
492!> \param atom_list list of selected atoms
493!> \param atom_map atomic map between device and contact force environments
494!> \param para_env parallel environment
495! **************************************************************************************************
496 SUBROUTINE negf_reference_contact_matrix(matrix_contact, matrix_device, atom_list, atom_map, para_env)
497 TYPE(dbcsr_type), POINTER :: matrix_contact, matrix_device
498 INTEGER, DIMENSION(:), INTENT(in) :: atom_list
499 TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
500 TYPE(mp_para_env_type), POINTER :: para_env
501
502 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_reference_contact_matrix'
503
504 INTEGER :: handle, i1, i2, iatom_col, iatom_row, &
505 icol, iproc, irow, max_atom, &
506 mepos_plus1, n1, n2, natoms, offset
507 INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_nelems, send_nelems
508 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: rank_contact, rank_device
509 LOGICAL :: found, transp
510 REAL(kind=dp), DIMENSION(:, :), POINTER :: rblock
511 TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_handlers, send_handlers
512 TYPE(negf_allocatable_rvector), ALLOCATABLE, &
513 DIMENSION(:) :: recv_packed_blocks, send_packed_blocks
514
515 CALL timeset(routinen, handle)
516 mepos_plus1 = para_env%mepos + 1
517
518 natoms = SIZE(atom_list)
519 max_atom = 0
520 DO iatom_row = 1, natoms
521 IF (atom_map(iatom_row)%iatom > max_atom) max_atom = atom_map(iatom_row)%iatom
522 END DO
523
524 ! find out which block goes to which node
525 ALLOCATE (rank_contact(max_atom, max_atom))
526 ALLOCATE (rank_device(max_atom, max_atom))
527
528 rank_contact(:, :) = 0
529 rank_device(:, :) = 0
530
531 DO iatom_col = 1, natoms
532 DO iatom_row = 1, iatom_col
533 IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
534 icol = atom_map(iatom_col)%iatom
535 irow = atom_map(iatom_row)%iatom
536 ELSE
537 icol = atom_map(iatom_row)%iatom
538 irow = atom_map(iatom_col)%iatom
539 END IF
540
541 CALL dbcsr_get_block_p(matrix=matrix_device, &
542 row=atom_list(iatom_row), col=atom_list(iatom_col), &
543 block=rblock, found=found)
544 IF (found) rank_device(irow, icol) = mepos_plus1
545
546 CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
547 IF (found) rank_contact(irow, icol) = mepos_plus1
548 END DO
549 END DO
550
551 CALL para_env%sum(rank_device)
552 CALL para_env%sum(rank_contact)
553
554 ! compute number of packed matrix elements to send to / receive from each processor
555 ALLOCATE (recv_nelems(para_env%num_pe))
556 ALLOCATE (send_nelems(para_env%num_pe))
557 recv_nelems(:) = 0
558 send_nelems(:) = 0
559
560 DO iatom_col = 1, natoms
561 DO iatom_row = 1, iatom_col
562 IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
563 icol = atom_map(iatom_col)%iatom
564 irow = atom_map(iatom_row)%iatom
565 ELSE
566 icol = atom_map(iatom_row)%iatom
567 irow = atom_map(iatom_col)%iatom
568 END IF
569
570 CALL dbcsr_get_block_p(matrix=matrix_device, &
571 row=atom_list(iatom_row), col=atom_list(iatom_col), &
572 block=rblock, found=found)
573 IF (found) THEN
574 iproc = rank_contact(irow, icol)
575 IF (iproc > 0) &
576 send_nelems(iproc) = send_nelems(iproc) + SIZE(rblock)
577 END IF
578
579 CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
580 IF (found) THEN
581 iproc = rank_device(irow, icol)
582 IF (iproc > 0) &
583 recv_nelems(iproc) = recv_nelems(iproc) + SIZE(rblock)
584 END IF
585 END DO
586 END DO
587
588 ! pack blocks
589 ALLOCATE (recv_packed_blocks(para_env%num_pe))
590 DO iproc = 1, para_env%num_pe
591 IF (iproc /= mepos_plus1 .AND. recv_nelems(iproc) > 0) &
592 ALLOCATE (recv_packed_blocks(iproc)%vector(recv_nelems(iproc)))
593 END DO
594
595 ALLOCATE (send_packed_blocks(para_env%num_pe))
596 DO iproc = 1, para_env%num_pe
597 IF (send_nelems(iproc) > 0) &
598 ALLOCATE (send_packed_blocks(iproc)%vector(send_nelems(iproc)))
599 END DO
600
601 send_nelems(:) = 0
602 DO iatom_col = 1, natoms
603 DO iatom_row = 1, iatom_col
604 IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
605 icol = atom_map(iatom_col)%iatom
606 irow = atom_map(iatom_row)%iatom
607 transp = .false.
608 ELSE
609 icol = atom_map(iatom_row)%iatom
610 irow = atom_map(iatom_col)%iatom
611 transp = .true.
612 END IF
613
614 iproc = rank_contact(irow, icol)
615 IF (iproc > 0) THEN
616 CALL dbcsr_get_block_p(matrix=matrix_device, &
617 row=atom_list(iatom_row), col=atom_list(iatom_col), &
618 block=rblock, found=found)
619 IF (found) THEN
620 offset = send_nelems(iproc)
621 n1 = SIZE(rblock, 1)
622 n2 = SIZE(rblock, 2)
623
624 IF (transp) THEN
625 DO i1 = 1, n1
626 DO i2 = 1, n2
627 send_packed_blocks(iproc)%vector(offset + i2) = rblock(i1, i2)
628 END DO
629 offset = offset + n2
630 END DO
631 ELSE
632 DO i2 = 1, n2
633 DO i1 = 1, n1
634 send_packed_blocks(iproc)%vector(offset + i1) = rblock(i1, i2)
635 END DO
636 offset = offset + n1
637 END DO
638 END IF
639
640 send_nelems(iproc) = offset
641 END IF
642 END IF
643 END DO
644 END DO
645
646 ! send blocks
647 ALLOCATE (recv_handlers(para_env%num_pe), send_handlers(para_env%num_pe))
648
649 DO iproc = 1, para_env%num_pe
650 IF (iproc /= mepos_plus1 .AND. send_nelems(iproc) > 0) THEN
651 CALL para_env%isend(send_packed_blocks(iproc)%vector, iproc - 1, send_handlers(iproc), 1)
652 END IF
653 END DO
654
655 ! receive blocks
656 DO iproc = 1, para_env%num_pe
657 IF (iproc /= mepos_plus1) THEN
658 IF (recv_nelems(iproc) > 0) THEN
659 CALL para_env%irecv(recv_packed_blocks(iproc)%vector, iproc - 1, recv_handlers(iproc), 1)
660 END IF
661 ELSE
662 IF (ALLOCATED(send_packed_blocks(iproc)%vector)) &
663 CALL move_alloc(send_packed_blocks(iproc)%vector, recv_packed_blocks(iproc)%vector)
664 END IF
665 END DO
666
667 ! unpack blocks
668 DO iproc = 1, para_env%num_pe
669 IF (iproc /= mepos_plus1 .AND. recv_nelems(iproc) > 0) &
670 CALL recv_handlers(iproc)%wait()
671 END DO
672
673 recv_nelems(:) = 0
674 DO iatom_col = 1, natoms
675 DO iatom_row = 1, iatom_col
676 IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
677 icol = atom_map(iatom_col)%iatom
678 irow = atom_map(iatom_row)%iatom
679 ELSE
680 icol = atom_map(iatom_row)%iatom
681 irow = atom_map(iatom_col)%iatom
682 END IF
683
684 iproc = rank_device(irow, icol)
685 IF (iproc > 0) THEN
686 CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
687
688 IF (found) THEN
689 offset = recv_nelems(iproc)
690 n1 = SIZE(rblock, 1)
691 n2 = SIZE(rblock, 2)
692
693 DO i2 = 1, n2
694 DO i1 = 1, n1
695 rblock(i1, i2) = recv_packed_blocks(iproc)%vector(offset + i1)
696 END DO
697 offset = offset + n1
698 END DO
699
700 recv_nelems(iproc) = offset
701 END IF
702 END IF
703 END DO
704 END DO
705
706 DO iproc = 1, para_env%num_pe
707 IF (iproc /= mepos_plus1 .AND. send_nelems(iproc) > 0) &
708 CALL send_handlers(iproc)%wait()
709 END DO
710
711 ! release memory
712 DEALLOCATE (recv_handlers, send_handlers)
713
714 DO iproc = para_env%num_pe, 1, -1
715 IF (ALLOCATED(send_packed_blocks(iproc)%vector)) &
716 DEALLOCATE (send_packed_blocks(iproc)%vector)
717 END DO
718 DEALLOCATE (send_packed_blocks)
719
720 DO iproc = para_env%num_pe, 1, -1
721 IF (ALLOCATED(recv_packed_blocks(iproc)%vector)) &
722 DEALLOCATE (recv_packed_blocks(iproc)%vector)
723 END DO
724 DEALLOCATE (recv_packed_blocks)
725
726 DEALLOCATE (rank_contact, rank_device)
727 CALL timestop(handle)
728 END SUBROUTINE negf_reference_contact_matrix
729
730! **************************************************************************************************
731!> \brief Invert cell_to_index mapping between unit cells and DBCSR matrix images.
732!> \param cell_to_index mapping: unit_cell -> image_index
733!> \param nimages number of images
734!> \param index_to_cell inverted mapping: image_index -> unit_cell
735!> \par History
736!> * 10.2017 created [Sergey Chulkov]
737! **************************************************************************************************
738 SUBROUTINE invert_cell_to_index(cell_to_index, nimages, index_to_cell)
739 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
740 INTEGER, INTENT(in) :: nimages
741 INTEGER, DIMENSION(3, nimages), INTENT(out) :: index_to_cell
742
743 CHARACTER(LEN=*), PARAMETER :: routinen = 'invert_cell_to_index'
744
745 INTEGER :: handle, i1, i2, i3, image
746 INTEGER, DIMENSION(3) :: lbounds, ubounds
747
748 CALL timeset(routinen, handle)
749
750 index_to_cell(:, :) = 0
751 lbounds = lbound(cell_to_index)
752 ubounds = ubound(cell_to_index)
753
754 DO i3 = lbounds(3), ubounds(3) ! z
755 DO i2 = lbounds(2), ubounds(2) ! y
756 DO i1 = lbounds(1), ubounds(1) ! x
757 image = cell_to_index(i1, i2, i3)
758 IF (image > 0 .AND. image <= nimages) THEN
759 index_to_cell(1, image) = i1
760 index_to_cell(2, image) = i2
761 index_to_cell(3, image) = i3
762 END IF
763 END DO
764 END DO
765 END DO
766
767 CALL timestop(handle)
768 END SUBROUTINE invert_cell_to_index
769
770! **************************************************************************************************
771!> \brief Helper routine to obtain index of a DBCSR matrix image by its unit cell replica.
772!> Can be used with any usin cell.
773!> \param cell indices of the unit cell
774!> \param cell_to_index mapping: unit_cell -> image_index
775!> \return DBCSR matrix images
776!> (0 means there are no non-zero matrix elements in the image)
777!> \par History
778!> * 10.2017 created [Sergey Chulkov]
779! **************************************************************************************************
780 PURE FUNCTION get_index_by_cell(cell, cell_to_index) RESULT(image)
781 INTEGER, DIMENSION(3), INTENT(in) :: cell
782 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
783 INTEGER :: image
784
785 IF (lbound(cell_to_index, 1) <= cell(1) .AND. ubound(cell_to_index, 1) >= cell(1) .AND. &
786 lbound(cell_to_index, 2) <= cell(2) .AND. ubound(cell_to_index, 2) >= cell(2) .AND. &
787 lbound(cell_to_index, 3) <= cell(3) .AND. ubound(cell_to_index, 3) >= cell(3)) THEN
788
789 image = cell_to_index(cell(1), cell(2), cell(3))
790 ELSE
791 image = 0
792 END IF
793 END FUNCTION get_index_by_cell
794END MODULE negf_matrix_utils
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:167
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1016
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:768
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:901
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
negf_alloc_types
Allocatable vectors for NEGF based quantum transport calculations.
Definition negf_alloc_types.F:12
negf_atom_map
Map atoms between various force environments.
Definition negf_atom_map.F:13
negf_matrix_utils
Helper routines to manipulate with matrices.
Definition negf_matrix_utils.F:12
negf_matrix_utils::negf_reference_contact_matrix
subroutine, public negf_reference_contact_matrix(matrix_contact, matrix_device, atom_list, atom_map, para_env)
Extract part of the DBCSR matrix based on selected atoms and copy it into another DBCSR matrix.
Definition negf_matrix_utils.F:497
negf_matrix_utils::number_of_atomic_orbitals
integer function, public number_of_atomic_orbitals(subsys, atom_list)
Compute the number of atomic orbitals of the given set of atoms.
Definition negf_matrix_utils.F:55
negf_matrix_utils::get_index_by_cell
pure integer function, public get_index_by_cell(cell, cell_to_index)
Helper routine to obtain index of a DBCSR matrix image by its unit cell replica. Can be used with any...
Definition negf_matrix_utils.F:781
negf_matrix_utils::invert_cell_to_index
subroutine, public invert_cell_to_index(cell_to_index, nimages, index_to_cell)
Invert cell_to_index mapping between unit cells and DBCSR matrix images.
Definition negf_matrix_utils.F:739
negf_matrix_utils::negf_copy_fm_submat_to_dbcsr
subroutine, public negf_copy_fm_submat_to_dbcsr(fm, matrix, atomlist_row, atomlist_col, subsys)
Populate relevant blocks of the DBCSR matrix using data from a ScaLAPACK matrix. Irrelevant blocks of...
Definition negf_matrix_utils.F:94
negf_matrix_utils::negf_copy_contact_matrix
subroutine, public negf_copy_contact_matrix(fm_cell0, fm_cell1, direction_axis, matrix_kp, index_to_cell, atom_list0, atom_list1, subsys, mpi_comm_global, is_same_cell, matrix_ref)
Driver routine to extract diagonal and off-diagonal blocks from a symmetric DBCSR matrix.
Definition negf_matrix_utils.F:297
negf_matrix_utils::negf_copy_sym_dbcsr_to_fm_submat
subroutine, public negf_copy_sym_dbcsr_to_fm_submat(matrix, fm, atomlist_row, atomlist_col, subsys, mpi_comm_global, do_upper_diag, do_lower)
Extract part of the DBCSR matrix based on selected atoms and copy it into a dense matrix.
Definition negf_matrix_utils.F:181
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition qs_subsys_types.F:134
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
message_passing::mp_request_type
Definition message_passing.F:623
negf_alloc_types::negf_allocatable_rvector
Allocatable 1-D real vector.
Definition negf_alloc_types.F:35
negf_atom_map::negf_atom_map_type
Structure that maps the given atom in the sourse FORCE_EVAL section with another atom from the target...
Definition negf_atom_map.F:38
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56