Line data Source code
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 :
12 : MODULE negf_matrix_utils
13 : USE cp_dbcsr_api, ONLY: &
14 : dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_get_block_p, dbcsr_get_info, &
15 : dbcsr_init_p, dbcsr_p_type, dbcsr_set, dbcsr_type
16 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
17 : USE cp_fm_types, ONLY: cp_fm_get_info,&
18 : cp_fm_get_submatrix,&
19 : cp_fm_set_submatrix,&
20 : cp_fm_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 :
44 : CONTAINS
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 24 : 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 24 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
62 24 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
63 :
64 24 : CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
65 72 : ALLOCATE (nsgfs(SIZE(particle_set)))
66 24 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
67 :
68 24 : IF (PRESENT(atom_list)) THEN
69 24 : natoms = SIZE(atom_list)
70 24 : nao = 0
71 :
72 152 : DO iatom = 1, natoms
73 152 : nao = nao + nsgfs(atom_list(iatom))
74 : END DO
75 : ELSE
76 0 : nao = SUM(nsgfs)
77 : END IF
78 :
79 24 : DEALLOCATE (nsgfs)
80 24 : 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 0 : 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 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
107 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
109 :
110 0 : CALL timeset(routineN, handle)
111 :
112 0 : CPASSERT(ASSOCIATED(matrix))
113 0 : CPASSERT(ASSOCIATED(subsys))
114 :
115 0 : CALL cp_fm_get_info(fm, nrow_global=nrows, ncol_global=ncols)
116 :
117 0 : CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
118 :
119 0 : natoms_row = SIZE(atomlist_row)
120 0 : natoms_col = SIZE(atomlist_col)
121 0 : nparticles = SIZE(particle_set)
122 :
123 0 : ALLOCATE (nsgfs(nparticles))
124 0 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
125 :
126 0 : ALLOCATE (fm_block(nrows, ncols))
127 0 : CALL cp_fm_get_submatrix(fm, fm_block)
128 :
129 0 : first_sgf_col = 1
130 0 : DO iatom_col = 1, natoms_col
131 : first_sgf_row = 1
132 0 : 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 0 : block=sm_block, found=found)
135 0 : 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 0 : nrows = nsgfs(atomlist_row(iatom_row))
140 0 : ncols = nsgfs(atomlist_col(iatom_col))
141 0 : DO icol = 1, ncols
142 0 : DO irow = 1, nrows
143 0 : 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 0 : first_sgf_row = first_sgf_row + nsgfs(atomlist_row(iatom_row))
149 : END DO
150 0 : first_sgf_col = first_sgf_col + nsgfs(atomlist_col(iatom_col))
151 : END DO
152 :
153 0 : DEALLOCATE (fm_block)
154 0 : DEALLOCATE (nsgfs)
155 :
156 0 : CALL timestop(handle)
157 0 : 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 96 : 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 96 : REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block
197 : TYPE(mp_para_env_type), POINTER :: para_env
198 96 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
199 96 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
200 :
201 96 : CALL timeset(routineN, handle)
202 :
203 96 : CPASSERT(ASSOCIATED(matrix))
204 96 : CPASSERT(ASSOCIATED(subsys))
205 :
206 96 : CALL qs_subsys_get(subsys, particle_set=particle_set, qs_kind_set=qs_kind_set)
207 :
208 96 : natoms_row = SIZE(atomlist_row)
209 96 : natoms_col = SIZE(atomlist_col)
210 96 : nparticles = SIZE(particle_set)
211 :
212 288 : ALLOCATE (nsgfs(nparticles))
213 96 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=nsgfs)
214 :
215 96 : 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 768 : CPASSERT(SUM(nsgfs(atomlist_row(:))) == nrows_fm)
219 640 : CPASSERT(SUM(nsgfs(atomlist_col(:))) == ncols_fm)
220 : END IF
221 :
222 384 : ALLOCATE (r2d(nrows_fm, ncols_fm))
223 19616 : r2d(:, :) = 0.0_dp
224 :
225 : offset_sgf_col = 0
226 640 : DO iatom_col = 1, natoms_col
227 : offset_sgf_row = 0
228 :
229 5152 : DO iatom_row = 1, natoms_row
230 4608 : IF (atomlist_row(iatom_row) <= atomlist_col(iatom_col)) THEN
231 2520 : IF (do_upper_diag) THEN
232 : CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_row(iatom_row), col=atomlist_col(iatom_col), &
233 2400 : block=sm_block, found=found)
234 : END IF
235 : ELSE
236 2088 : IF (do_lower) THEN
237 : CALL dbcsr_get_block_p(matrix=matrix, row=atomlist_col(iatom_col), col=atomlist_row(iatom_row), &
238 2016 : block=sm_block, found=found)
239 : END IF
240 : END IF
241 :
242 4608 : IF (found) THEN
243 1588 : IF (atomlist_row(iatom_row) <= atomlist_col(iatom_col)) THEN
244 899 : IF (do_upper_diag) THEN
245 2535 : DO icol = nsgfs(atomlist_col(iatom_col)), 1, -1
246 5915 : DO irow = nsgfs(atomlist_row(iatom_row)), 1, -1
247 5070 : 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 689 : IF (do_lower) THEN
253 1959 : DO icol = nsgfs(atomlist_col(iatom_col)), 1, -1
254 4571 : DO irow = nsgfs(atomlist_row(iatom_row)), 1, -1
255 3918 : 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 5152 : offset_sgf_row = offset_sgf_row + nsgfs(atomlist_row(iatom_row))
263 : END DO
264 640 : offset_sgf_col = offset_sgf_col + nsgfs(atomlist_col(iatom_col))
265 : END DO
266 :
267 96 : CALL mpi_comm_global%sum(r2d)
268 :
269 96 : CALL cp_fm_set_submatrix(fm, r2d)
270 :
271 96 : DEALLOCATE (r2d)
272 96 : DEALLOCATE (nsgfs)
273 :
274 96 : CALL timestop(handle)
275 192 : 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 12 : SUBROUTINE negf_copy_contact_matrix(fm_cell0, fm_cell1, direction_axis, matrix_kp, index_to_cell, &
296 12 : 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 12 : REAL(kind=dp), DIMENSION(:, :), POINTER :: block_dest, block_src
316 12 : 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 12 : CALL timeset(routineN, handle)
321 :
322 12 : CPASSERT(ASSOCIATED(subsys))
323 :
324 12 : nimages = SIZE(index_to_cell, 2)
325 12 : 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 84 : ALLOCATE (matrix_cells_raw(-2:2))
333 72 : DO rep = -2, 2
334 60 : NULLIFY (matrix_cells_raw(rep)%matrix)
335 60 : CALL dbcsr_init_p(matrix_cells_raw(rep)%matrix)
336 60 : CALL dbcsr_copy(matrix_cells_raw(rep)%matrix, matrix_kp(1)%matrix)
337 72 : CALL dbcsr_set(matrix_cells_raw(rep)%matrix, 0.0_dp)
338 : END DO
339 :
340 12 : NULLIFY (matrix_cell_0, matrix_cell_1, matrix_cell_minus1)
341 :
342 12 : CALL dbcsr_init_p(matrix_cell_0)
343 12 : CALL dbcsr_copy(matrix_cell_0, matrix_kp(1)%matrix)
344 12 : CALL dbcsr_set(matrix_cell_0, 0.0_dp)
345 :
346 12 : CALL dbcsr_init_p(matrix_cell_1)
347 12 : CALL dbcsr_copy(matrix_cell_1, matrix_kp(1)%matrix)
348 12 : CALL dbcsr_set(matrix_cell_1, 0.0_dp)
349 :
350 12 : CALL dbcsr_init_p(matrix_cell_minus1)
351 12 : CALL dbcsr_copy(matrix_cell_minus1, matrix_kp(1)%matrix)
352 12 : CALL dbcsr_set(matrix_cell_minus1, 0.0_dp)
353 :
354 1164 : DO image = 1, nimages
355 1152 : rep = index_to_cell(direction_axis_abs, image)
356 :
357 1152 : IF (ABS(rep) <= 2) &
358 1164 : CALL dbcsr_add(matrix_cells_raw(rep)%matrix, matrix_kp(image)%matrix, 1.0_dp, 1.0_dp)
359 : END DO
360 :
361 12 : CALL dbcsr_get_info(matrix_cell_0, nblkrows_total=natoms)
362 :
363 12 : 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 84 : is_same_cell(:, :) = 0
367 :
368 20 : DO iatom_col = 1, natoms
369 60 : DO iatom_row = 1, iatom_col
370 : CALL dbcsr_get_block_p(matrix=matrix_ref, &
371 : row=iatom_row, col=iatom_col, &
372 40 : block=block_src, found=found)
373 :
374 56 : 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 20 : IF (MOD(iatom_col - iatom_row, 2) == 0) THEN
378 : phase = 1
379 : ELSE
380 8 : 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 20 : block=block_dest, found=found)
386 20 : CPASSERT(found)
387 :
388 140 : 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 20 : block=block_dest, found=found)
393 20 : CPASSERT(found)
394 140 : error_diff = MAXVAL(ABS(block_dest(:, :) - block_src(:, :)))
395 :
396 20 : IF (error_same <= error_diff) THEN
397 14 : is_same_cell(iatom_row, iatom_col) = 0
398 : ELSE
399 6 : 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 60 : DO iatom_col = 1, natoms
407 180 : DO iatom_row = 1, iatom_col
408 : CALL dbcsr_get_block_p(matrix=matrix_cell_0, &
409 120 : row=iatom_row, col=iatom_col, block=block_dest, found=found)
410 :
411 168 : IF (found) THEN
412 : ! it should be much safe to rely on a neighbour list iterator
413 60 : IF (MOD(iatom_col - iatom_row, 2) == 0) THEN
414 : phase = 1
415 : ELSE
416 24 : phase = -1
417 : END IF
418 60 : 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 60 : row=iatom_row, col=iatom_col, block=block_src, found=found)
425 60 : CPASSERT(found)
426 840 : 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 60 : row=iatom_row, col=iatom_col, block=block_dest, found=found)
433 60 : CPASSERT(found)
434 : CALL dbcsr_get_block_p(matrix=matrix_cells_raw(rep + phase)%matrix, &
435 60 : row=iatom_row, col=iatom_col, block=block_src, found=found)
436 60 : CPASSERT(found)
437 840 : 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 60 : row=iatom_row, col=iatom_col, block=block_dest, found=found)
444 60 : CPASSERT(found)
445 : CALL dbcsr_get_block_p(matrix=matrix_cells_raw(rep - phase)%matrix, &
446 60 : row=iatom_row, col=iatom_col, block=block_src, found=found)
447 60 : CPASSERT(found)
448 840 : block_dest(:, :) = block_src(:, :)
449 : END IF
450 : END DO
451 : END DO
452 :
453 12 : 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 6 : 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 6 : 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 6 : 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 6 : subsys, mpi_comm_global, do_upper_diag=.FALSE., do_lower=.TRUE.)
467 :
468 : END IF
469 12 : 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 12 : subsys, mpi_comm_global, do_upper_diag=.TRUE., do_lower=.TRUE.)
474 :
475 12 : CALL dbcsr_deallocate_matrix(matrix_cell_0)
476 12 : CALL dbcsr_deallocate_matrix(matrix_cell_1)
477 12 : CALL dbcsr_deallocate_matrix(matrix_cell_minus1)
478 :
479 72 : DO rep = -2, 2
480 72 : CALL dbcsr_deallocate_matrix(matrix_cells_raw(rep)%matrix)
481 : END DO
482 12 : DEALLOCATE (matrix_cells_raw)
483 :
484 12 : CALL timestop(handle)
485 24 : 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 4 : 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 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: recv_nelems, send_nelems
508 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: rank_contact, rank_device
509 : LOGICAL :: found, transp
510 4 : REAL(kind=dp), DIMENSION(:, :), POINTER :: rblock
511 4 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_handlers, send_handlers
512 : TYPE(negf_allocatable_rvector), ALLOCATABLE, &
513 4 : DIMENSION(:) :: recv_packed_blocks, send_packed_blocks
514 :
515 4 : CALL timeset(routineN, handle)
516 4 : mepos_plus1 = para_env%mepos + 1
517 :
518 4 : natoms = SIZE(atom_list)
519 4 : max_atom = 0
520 20 : DO iatom_row = 1, natoms
521 20 : 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 16 : ALLOCATE (rank_contact(max_atom, max_atom))
526 12 : ALLOCATE (rank_device(max_atom, max_atom))
527 :
528 84 : rank_contact(:, :) = 0
529 84 : rank_device(:, :) = 0
530 :
531 20 : DO iatom_col = 1, natoms
532 60 : DO iatom_row = 1, iatom_col
533 40 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
534 40 : icol = atom_map(iatom_col)%iatom
535 40 : irow = atom_map(iatom_row)%iatom
536 : ELSE
537 0 : icol = atom_map(iatom_row)%iatom
538 0 : 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 40 : block=rblock, found=found)
544 40 : IF (found) rank_device(irow, icol) = mepos_plus1
545 :
546 40 : CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
547 96 : IF (found) rank_contact(irow, icol) = mepos_plus1
548 : END DO
549 : END DO
550 :
551 4 : CALL para_env%sum(rank_device)
552 4 : CALL para_env%sum(rank_contact)
553 :
554 : ! compute number of packed matrix elements to send to / receive from each processor
555 12 : ALLOCATE (recv_nelems(para_env%num_pe))
556 8 : ALLOCATE (send_nelems(para_env%num_pe))
557 12 : recv_nelems(:) = 0
558 12 : send_nelems(:) = 0
559 :
560 20 : DO iatom_col = 1, natoms
561 60 : DO iatom_row = 1, iatom_col
562 40 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
563 40 : icol = atom_map(iatom_col)%iatom
564 40 : irow = atom_map(iatom_row)%iatom
565 : ELSE
566 0 : icol = atom_map(iatom_row)%iatom
567 0 : 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 40 : block=rblock, found=found)
573 40 : IF (found) THEN
574 20 : iproc = rank_contact(irow, icol)
575 20 : IF (iproc > 0) &
576 60 : send_nelems(iproc) = send_nelems(iproc) + SIZE(rblock)
577 : END IF
578 :
579 40 : CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
580 96 : IF (found) THEN
581 20 : iproc = rank_device(irow, icol)
582 20 : IF (iproc > 0) &
583 60 : recv_nelems(iproc) = recv_nelems(iproc) + SIZE(rblock)
584 : END IF
585 : END DO
586 : END DO
587 :
588 : ! pack blocks
589 20 : ALLOCATE (recv_packed_blocks(para_env%num_pe))
590 12 : DO iproc = 1, para_env%num_pe
591 8 : IF (iproc /= mepos_plus1 .AND. recv_nelems(iproc) > 0) &
592 4 : ALLOCATE (recv_packed_blocks(iproc)%vector(recv_nelems(iproc)))
593 : END DO
594 :
595 20 : ALLOCATE (send_packed_blocks(para_env%num_pe))
596 12 : DO iproc = 1, para_env%num_pe
597 8 : IF (send_nelems(iproc) > 0) &
598 16 : ALLOCATE (send_packed_blocks(iproc)%vector(send_nelems(iproc)))
599 : END DO
600 :
601 12 : send_nelems(:) = 0
602 20 : DO iatom_col = 1, natoms
603 60 : DO iatom_row = 1, iatom_col
604 40 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
605 40 : icol = atom_map(iatom_col)%iatom
606 40 : irow = atom_map(iatom_row)%iatom
607 40 : transp = .FALSE.
608 : ELSE
609 0 : icol = atom_map(iatom_row)%iatom
610 0 : irow = atom_map(iatom_col)%iatom
611 0 : transp = .TRUE.
612 : END IF
613 :
614 40 : iproc = rank_contact(irow, icol)
615 56 : 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 40 : block=rblock, found=found)
619 40 : IF (found) THEN
620 20 : offset = send_nelems(iproc)
621 20 : n1 = SIZE(rblock, 1)
622 20 : n2 = SIZE(rblock, 2)
623 :
624 20 : IF (transp) THEN
625 0 : DO i1 = 1, n1
626 0 : DO i2 = 1, n2
627 0 : send_packed_blocks(iproc)%vector(offset + i2) = rblock(i1, i2)
628 : END DO
629 0 : offset = offset + n2
630 : END DO
631 : ELSE
632 60 : DO i2 = 1, n2
633 120 : DO i1 = 1, n1
634 120 : send_packed_blocks(iproc)%vector(offset + i1) = rblock(i1, i2)
635 : END DO
636 60 : offset = offset + n1
637 : END DO
638 : END IF
639 :
640 20 : send_nelems(iproc) = offset
641 : END IF
642 : END IF
643 : END DO
644 : END DO
645 :
646 : ! send blocks
647 36 : ALLOCATE (recv_handlers(para_env%num_pe), send_handlers(para_env%num_pe))
648 :
649 12 : DO iproc = 1, para_env%num_pe
650 12 : IF (iproc /= mepos_plus1 .AND. send_nelems(iproc) > 0) THEN
651 0 : 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 12 : DO iproc = 1, para_env%num_pe
657 12 : IF (iproc /= mepos_plus1) THEN
658 4 : IF (recv_nelems(iproc) > 0) THEN
659 0 : CALL para_env%irecv(recv_packed_blocks(iproc)%vector, iproc - 1, recv_handlers(iproc), 1)
660 : END IF
661 : ELSE
662 4 : IF (ALLOCATED(send_packed_blocks(iproc)%vector)) &
663 4 : CALL MOVE_ALLOC(send_packed_blocks(iproc)%vector, recv_packed_blocks(iproc)%vector)
664 : END IF
665 : END DO
666 :
667 : ! unpack blocks
668 12 : DO iproc = 1, para_env%num_pe
669 8 : IF (iproc /= mepos_plus1 .AND. recv_nelems(iproc) > 0) &
670 4 : CALL recv_handlers(iproc)%wait()
671 : END DO
672 :
673 12 : recv_nelems(:) = 0
674 20 : DO iatom_col = 1, natoms
675 60 : DO iatom_row = 1, iatom_col
676 40 : IF (atom_map(iatom_row)%iatom <= atom_map(iatom_col)%iatom) THEN
677 40 : icol = atom_map(iatom_col)%iatom
678 40 : irow = atom_map(iatom_row)%iatom
679 : ELSE
680 0 : icol = atom_map(iatom_row)%iatom
681 0 : irow = atom_map(iatom_col)%iatom
682 : END IF
683 :
684 40 : iproc = rank_device(irow, icol)
685 56 : IF (iproc > 0) THEN
686 40 : CALL dbcsr_get_block_p(matrix=matrix_contact, row=irow, col=icol, block=rblock, found=found)
687 :
688 40 : IF (found) THEN
689 20 : offset = recv_nelems(iproc)
690 20 : n1 = SIZE(rblock, 1)
691 20 : n2 = SIZE(rblock, 2)
692 :
693 60 : DO i2 = 1, n2
694 120 : DO i1 = 1, n1
695 120 : rblock(i1, i2) = recv_packed_blocks(iproc)%vector(offset + i1)
696 : END DO
697 60 : offset = offset + n1
698 : END DO
699 :
700 20 : recv_nelems(iproc) = offset
701 : END IF
702 : END IF
703 : END DO
704 : END DO
705 :
706 12 : DO iproc = 1, para_env%num_pe
707 8 : IF (iproc /= mepos_plus1 .AND. send_nelems(iproc) > 0) &
708 4 : CALL send_handlers(iproc)%wait()
709 : END DO
710 :
711 : ! release memory
712 4 : DEALLOCATE (recv_handlers, send_handlers)
713 :
714 12 : DO iproc = para_env%num_pe, 1, -1
715 8 : IF (ALLOCATED(send_packed_blocks(iproc)%vector)) &
716 4 : DEALLOCATE (send_packed_blocks(iproc)%vector)
717 : END DO
718 12 : DEALLOCATE (send_packed_blocks)
719 :
720 12 : DO iproc = para_env%num_pe, 1, -1
721 8 : IF (ALLOCATED(recv_packed_blocks(iproc)%vector)) &
722 8 : DEALLOCATE (recv_packed_blocks(iproc)%vector)
723 : END DO
724 12 : DEALLOCATE (recv_packed_blocks)
725 :
726 4 : DEALLOCATE (rank_contact, rank_device)
727 4 : CALL timestop(handle)
728 8 : 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 8 : 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 8 : CALL timeset(routineN, handle)
749 :
750 3080 : index_to_cell(:, :) = 0
751 32 : lbounds = LBOUND(cell_to_index)
752 32 : ubounds = UBOUND(cell_to_index)
753 :
754 48 : DO i3 = lbounds(3), ubounds(3) ! z
755 328 : DO i2 = lbounds(2), ubounds(2) ! y
756 2280 : DO i1 = lbounds(1), ubounds(1) ! x
757 1960 : image = cell_to_index(i1, i2, i3)
758 2240 : IF (image > 0 .AND. image <= nimages) THEN
759 768 : index_to_cell(1, image) = i1
760 768 : index_to_cell(2, image) = i2
761 768 : index_to_cell(3, image) = i3
762 : END IF
763 : END DO
764 : END DO
765 : END DO
766 :
767 8 : CALL timestop(handle)
768 8 : 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 7040 : 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 49280 : LBOUND(cell_to_index, 3) <= cell(3) .AND. UBOUND(cell_to_index, 3) >= cell(3)) THEN
788 :
789 7040 : image = cell_to_index(cell(1), cell(2), cell(3))
790 : ELSE
791 : image = 0
792 : END IF
793 7040 : END FUNCTION get_index_by_cell
794 : END MODULE negf_matrix_utils
|
