1!--------------------------------------------------------------------------------------------------! 2! DFTB+: general package for performing fast atomistic simulations ! 3! Copyright (C) 2006 - 2019 DFTB+ developers group ! 4! ! 5! See the LICENSE file for terms of usage and distribution. ! 6!--------------------------------------------------------------------------------------------------! 7 8 9#:include "common.fypp" 10 11 12!> Contains routines for folding internal sparse matrixes into Compressed 13!> Sparse Row Format and vica-versa. 14!> 15!> Caveat All folding and unfolding routines assume, that the CSR matrices have the same sparsity 16!> structure as the internal sparse matrices. The CSR matrices must be created using the 17!> createEquivCSR routines. 18module dftbp_matconv 19 use dftbp_assert 20 use dftbp_accuracy 21 use dftbp_constants, only : pi 22 use dftbp_commontypes, only : TOrbitals 23 use libnegf, only: r_CSR, z_CSR, r_DNS, z_DNS, create, destroy 24 implicit none 25 private 26 27 public :: init, destruct 28 public :: foldToCSR, unfoldFromCSR 29 30 31 !> Creates a Compressed Sparse Row matrix with equivalent sparsity pattern to the internal sparse 32 !> matrices. 33 interface init 34 module procedure createEmptyCSR_real 35 module procedure createEmptyCSR_cplx 36 module procedure createEquivCSR_real 37 module procedure createEquivCSR_cplx 38 module procedure cloneSparsityMap_real_real 39 module procedure cloneSparsityMap_cplx_cplx 40 module procedure cloneSparsityMap_real_cplx 41 module procedure cloneSparsityMap_cplx_real 42 end interface init 43 44 !> Destructor for the different structures. 45 interface destruct 46 module procedure rdestroy_CSR, zdestroy_CSR 47 module procedure rdestroy_DNS, zdestroy_DNS 48 end interface destruct 49 50 !> Folds internal sparse storage to Compressed Sparse Row format. 51 interface foldToCSR 52 module procedure foldToCSR_real 53 module procedure foldToCSR_cplx 54 end interface foldToCSR 55 56 !> Unfolds from Compressed Sparse Row format into internal sparse storage. 57 interface unfoldFromCSR 58 module procedure unfoldFromCSR_real 59 module procedure unfoldFromCSR_cplx 60 end interface unfoldFromCSR 61 62contains 63 64 !> Creates a Compressed Sparse Row matrix with a sparsity structure in 65 !> accordance with current geometry. 66 !> 67 !> Note: The subroutine only creates the structure of the CSR matrix, but 68 !> does not fill up the nonzero values with anything usefull. 69 !> The resulting CSR matrix has storage for both triangles of the square matrix. 70 subroutine createEquivCSR_real(csr, iAtomStart, iNeighbor, nNeighbor, img2CentCell, orb) 71 72 !> Compressed sparse row matrix on exit 73 type(r_CSR), intent(out) :: csr 74 75 !> Starting position of the atoms in the square H/S. 76 integer, intent(in) :: iAtomStart(:) 77 78 !> Neighborlist for each atom. 79 integer, intent(in) :: iNeighbor(0:,:) 80 81 !> Number of neighbors for each atom. 82 integer, intent(in) :: nNeighbor(:) 83 84 !> Folded image in the central cell for each atom. 85 integer, intent(in) :: img2CentCell(:) 86 87 !> orb Orbitals in the system. 88 type(TOrbitals), intent(in) :: orb 89 90 integer, allocatable :: nColAtom(:), nCols(:), iNonZero(:) 91 integer, allocatable :: rowpnt(:) 92 logical, allocatable :: zero(:) 93 integer :: nAtom, nNonZero 94 integer :: iOrb1, iOrb2, nOrb1, nOrb2, iAt1, iAt2f, iNeigh 95 integer :: iRow, iCol, ind, nrow 96 integer :: ii, jj 97 98 nAtom = size(orb%nOrbAtom) 99 100 ! Count nr. of nonzero columns in the square (folded) form for each atom. Holds nr. of nonzero 101 ! columns for each atom. 102 allocate(nColAtom(nAtom)) 103 allocate(zero(nAtom)) 104 nColAtom(:) = 0 105 do iAt1 = 1, nAtom 106 zero(:) = .true. 107 nOrb1 = orb%nOrbAtom(iAt1) 108 do iNeigh = 0, nNeighbor(iAt1) 109 iAt2f = img2CentCell(iNeighbor(iNeigh,iAt1)) 110 if (zero(iAt2f)) then 111 nColAtom(iAt1) = nColAtom(iAt1) + orb%nOrbAtom(iAt2f) 112 zero(iAt2f) = .false. 113 if (iAt1 /= iAt2f) then 114 nColAtom(iAt2f) = nColAtom(iAt2f) + nOrb1 115 end if 116 end if 117 end do 118 end do 119 120 ! Calculate CSR row pointers: 121 ! A row for a certain orbital of a certain atom is as long, as the 122 ! nr. of columns determined previously for the atom. 123 nrow = iAtomStart(nAtom+1) - 1 124 allocate(rowpnt(nrow+1)) 125 rowpnt(1) = 1 126 ind = 1 127 do iAt1 = 1, nAtom 128 nOrb1 = orb%nOrbAtom(iAt1) 129 do iOrb1 = 1, nOrb1 130 rowpnt(ind+iOrb1) = rowpnt(ind) + iOrb1 * nColAtom(iAt1) 131 end do 132 ind = ind + nOrb1 133 end do 134 deallocate(nColAtom) 135 136 call create(csr, nrow, nrow, rowpnt(nrow + 1) - 1) 137 csr%rowpnt = rowpnt 138 deallocate(rowpnt) 139 140 ! Nr. of CSR columns already filled 141 allocate(nCols(csr%nRow)) 142 nCols(:) = 0 143 ! Index of the nonzero blocks 144 allocate(iNonZero(nAtom)) 145 146 ! Loop over all atoms (over all atomic columns in the rectangular picture) 147 lpAt1: do iAt1 = 1, nAtom 148 iCol = iAtomStart(iAt1) 149 ! Width of the atomic column 150 nOrb1 = orb%nOrbAtom(iAt1) 151 152 ! Mark nonzero blocks in the folded matrix for the current atom 153 nNonZero = 0 154 zero(:) = .true. 155 do iNeigh = 0, nNeighbor(iAt1) 156 iAt2f = img2CentCell(iNeighbor(iNeigh,iAt1)) 157 if (zero(iAt2f)) then 158 zero(iAt2f) = .false. 159 nNonZero = nNonZero + 1 160 iNonzero(nNonzero) = iAt2f 161 end if 162 end do 163 164 ! Initialise CSR internal pointers according the non-zero blocks 165 lpNonZero: do ii = 1, nNonZero 166 iAt2f = iNonZero(ii) 167 nOrb2 = orb%nOrbAtom(iAt2f) 168 iRow = iAtomStart(iAt2f) 169 170 ! Correspond rows of the current atomic block column to the appropriate 171 ! partial rows in the lower triangle of the CSR matrix. 172 do iOrb2 = 0, nOrb2 - 1 173 jj = iRow + iOrb2 174 ind = csr%rowpnt(jj) + nCols(jj) 175 do iOrb1 = 0, nOrb1 - 1 176 csr%colind(ind+iOrb1) = iCol + iOrb1 177 end do 178 nCols(jj) = nCols(jj) + nOrb1 179 end do 180 181 ! Correspond the columns of the current block to appropriate 182 ! partial rows in the upper triangle of the CSR matrix. 183 if (iAt2f /= iAt1) then 184 do iOrb1 = 0, nOrb1 - 1 185 jj = iCol + iOrb1 186 ind = csr%rowpnt(jj) + nCols(jj) 187 do iOrb2 = 0, nOrb2 - 1 188 csr%colind(ind+iOrb2) = iRow + iOrb2 189 end do 190 nCols(jj) = nCols(jj) + nOrb2 191 end do 192 end if 193 end do lpNonZero 194 end do lpAt1 195 196 deallocate(zero) 197 deallocate(nCols) 198 deallocate(iNonZero) 199 200 end subroutine createEquivCSR_real 201 202 203 204 !> Folds the internal sparse formatted matrix to the compressed sparse row 205 !> format (real version). 206 !> 207 !> Note: In the resulting CSR format both triangles of the matrix are set. 208 !> Caveat: The routine should only applied on CSR matrixes created by the 209 !> createEquiv_real subroutine, since it exploits the ordering of the 210 !> column indexes. 211 subroutine foldToCSR_real(csr, sparse, iAtomStart, iPair, iNeighbor, nNeighbor, img2CentCell, orb) 212 213 !> csr Resulting CSR matrix 214 type(r_CSR), intent(inout) :: csr 215 216 !> sparse The internal sparse matrix to fold. 217 real(dp), intent(in) :: sparse(:) 218 219 !> iAtomStart Starting positions of the atoms in the square matrix 220 integer, intent(in) :: iAtomStart(:) 221 222 !> iPair Starting position of atom-neighbor interaction in the sparse matrix. 223 integer, intent(in) :: iPair(0:,:) 224 225 !> iNeighbor Index of neighbors 226 integer, intent(in) :: iNeighbor(0:,:) 227 228 !> nNeighbor Number of neighbors 229 integer, intent(in) :: nNeighbor(:) 230 231 !> img2CentCell Image of the atoms in the central cell. 232 integer, intent(in) :: img2CentCell(:) 233 234 !> orb Orbitals in the system. 235 type(TOrbitals), intent(in) :: orb 236 237 integer, allocatable :: nCols(:) 238 real(dp), allocatable :: tmpCol(:,:) 239 integer :: nAtom 240 integer :: iOrb1, nOrb1, nOrb2, iAt1, iAt2f, iNeigh 241 integer :: iRow, iCol, iVal, ind 242 integer :: ii, jj, kk 243 244 nAtom = size(orb%nOrbAtom) 245 246 csr%sorted = .false. 247 248 ! One atomic block column. 249 allocate(tmpCol(csr%nRow, orb%mOrb)) 250 lpAt1: do iAt1 = 1, nAtom 251 ! Unfold the columns for the current atom from the sparse matrix. 252 nOrb1 = orb%nOrbAtom(iAt1) 253 tmpCol(:,:) = 0.0_dp 254 do iNeigh = 0, nNeighbor(iAt1) 255 ind = iPair(iNeigh,iAt1) + 1 256 iAt2f = img2CentCell(iNeighbor(iNeigh,iAt1)) 257 nOrb2 = orb%nOrbAtom(iAt2f) 258 iRow = iAtomStart(iAt2f) 259 tmpCol(iRow:iRow+nOrb2-1, 1:nOrb1) = tmpCol(iRow:iRow+nOrb2-1, 1:nOrb1)& 260 & + reshape(sparse(ind:ind+nOrb2*nOrb1-1), (/ nOrb2, nOrb1 /)) 261 end do 262 263 ! Copy every column into the appropriate row in the upper triangle of 264 ! the CSR matrix (the copy is masked by the sparsity structure stored 265 ! in the CSR matrix) 266 nOrb1 = orb%nOrbAtom(iAt1) 267 iCol = iAtomStart(iAt1) 268 do iOrb1 = 1, nOrb1 269 ii = csr%rowpnt(iCol + iOrb1 - 1) 270 jj = csr%rowpnt(iCol + iOrb1) - 1 271 do kk=ii,jj 272 csr%nzval(kk) = tmpCol(csr%colind(kk), iOrb1) 273 enddo 274 end do 275 end do lpAt1 276 deallocate(tmpCol) 277 278 ! Fill the lower triangle of the CSR matrix 279 allocate(nCols(csr%nRow)) 280 nCols(:) = 0 281 do iRow = 1, csr%nRow 282 ! Starting from the tail of the matrix row 283 iVal = csr%rowpnt(iRow + 1) - 1 284 iCol = csr%colind(iVal) 285 do while (iVal >= csr%rowpnt(iRow) .and. iCol > iRow) 286 csr%nzval(csr%rowpnt(iCol)+nCols(iCol)) = csr%nzval(iVal) 287 nCols(iCol) = nCols(iCol) + 1 288 iVal = iVal - 1 289 iCol = csr%colind(iVal) 290 end do 291 end do 292 deallocate(nCols) 293 294 end subroutine foldToCSR_real 295 296 297 298 !> Unfolds a matrix from the CSR form into the internal DFTB+ sparse representation (real 299 !> version). 300 !> 301 !> Note: The CSR matrix must be symmetric. The unfolded matrix is added to the passed sparse 302 !> matrix. 303 subroutine unfoldFromCSR_real(sparse, csr, iAtomStart, iPair, iNeighbor, nNeighbor, img2CentCell,& 304 & orb) 305 306 !> sparse Updated sparse matrix. 307 real(dp), intent(inout) :: sparse(:) 308 309 !> CSR matrix 310 type(r_CSR), intent(in) :: csr 311 312 !> Starting positions of the atoms in the square matrix 313 integer, intent(in) :: iAtomStart(:) 314 315 !> Starting position of atom-neighbor interaction in the sparse matrix. 316 integer, intent(in) :: iPair(0:,:) 317 318 !> Index of neighbors 319 integer, intent(in) :: iNeighbor(0:,:) 320 321 !> Number of neighbors 322 integer, intent(in) :: nNeighbor(:) 323 324 !> Image of the atoms in the central cell. 325 integer, intent(in) :: img2CentCell(:) 326 327 !> Orbitals in the system. 328 type(TOrbitals), intent(in) :: orb 329 330 real(dp), allocatable :: tmpCol(:,:) 331 integer :: nAtom 332 integer :: iOrb1, nOrb1, nOrb2, iAt1, iAt2, iAt2f, iNeigh 333 integer :: iRow, iCol, ind 334 integer :: ii 335 336 nAtom = size(orb%nOrbAtom) 337 338 @:ASSERT(csr%nRow == iAtomStart(nAtom+1) - 1) 339 340 allocate(tmpCol(csr%nRow, orb%mOrb)) 341 342 do iAt1 = 1, nAtom 343 ! Put the rows belonging to a certain atom into the appropriate column 344 ! of the block column. (Matrix is assumed to be symmetric!) 345 tmpCol(:,:) = 0.0_dp 346 nOrb1 = orb%nOrbAtom(iAt1) 347 iCol = iAtomStart(iAt1) 348 do iOrb1 = 1, nOrb1 349 ii = iCol + iOrb1 - 1 350 do ind = csr%rowpnt(ii), csr%rowpnt(ii+1) - 1 351 tmpCol(csr%colind(ind), iOrb1) = real(csr%nzval(ind)) 352 end do 353 end do 354 355 ! Unfold the block column into the internal sparse format 356 do iNeigh = 0, nNeighbor(iAt1) 357 ind = iPair(iNeigh,iAt1) + 1 358 iAt2 = iNeighbor(iNeigh, iAt1) 359 iAt2f = img2CentCell(iAt2) 360 nOrb2 = orb%nOrbAtom(iAt2f) 361 iRow = iAtomStart(iAt2f) 362 sparse(ind:ind+nOrb2*nOrb1-1) = sparse(ind:ind+nOrb2*nOrb1-1)& 363 & + reshape(tmpCol(iRow:iRow+nOrb2-1, 1:nOrb1), (/nOrb2*nOrb1/)) 364 end do 365 end do 366 367 deallocate(tmpCol) 368 369 end subroutine unfoldFromCSR_real 370 371 372 373 !> Creates a Compressed Sparse Row matrix with a sparsity structure in accordance with current 374 !> geometry. 375 !> 376 !> Note: The subroutine only creates the structure of the CSR matrix, but does not fill up the 377 !> nonzero values with anything usefull. The resulting CSR matrix has storage for both triangles 378 !> of the square matrix. 379 subroutine createEquivCSR_cplx(csr, iAtomStart, iNeighbor, nNeighbor, img2CentCell, orb) 380 381 !> Compressed sparse row matrix on exit 382 type(z_CSR), intent(out) :: csr 383 384 !> Starting position of the atoms in the square H/S. 385 integer, intent(in) :: iAtomStart(:) 386 387 !> Neighborlist for each atom. 388 integer, intent(in) :: iNeighbor(0:,:) 389 390 !> Number of neighbors for each atom. 391 integer, intent(in) :: nNeighbor(:) 392 393 !> Folded image in the central cell for each atom. 394 integer, intent(in) :: img2CentCell(:) 395 396 !> Orbitals in the system. 397 type(TOrbitals), intent(in) :: orb 398 399 integer, allocatable :: nColAtom(:), nCols(:), iNonZero(:) 400 integer, allocatable :: rowpnt(:) 401 logical, allocatable :: zero(:) 402 integer :: nAtom, nNonZero 403 integer :: iOrb1, iOrb2, nOrb1, nOrb2, iAt1, iAt2f, iNeigh 404 integer :: iRow, iCol, ind, nrow 405 integer :: ii, jj 406 407 nAtom = size(orb%nOrbAtom) 408 409 ! Count nr. of nonzero columns in the square (folded) form for each atom 410 ! Nr. of nonzero columns for each atom 411 allocate(nColAtom(nAtom)) 412 allocate(zero(nAtom)) 413 nColAtom(:) = 0 414 do iAt1 = 1, nAtom 415 zero(:) = .true. 416 nOrb1 = orb%nOrbAtom(iAt1) 417 do iNeigh = 0, nNeighbor(iAt1) 418 iAt2f = img2CentCell(iNeighbor(iNeigh,iAt1)) 419 if (zero(iAt2f)) then 420 nColAtom(iAt1) = nColAtom(iAt1) + orb%nOrbAtom(iAt2f) 421 zero(iAt2f) = .false. 422 if (iAt1 /= iAt2f) then 423 nColAtom(iAt2f) = nColAtom(iAt2f) + nOrb1 424 end if 425 end if 426 end do 427 end do 428 429 nrow = iAtomStart(nAtom+1) - 1 430 allocate(rowpnt(nrow+1) ) 431 ! Calculate CSR row pointers: 432 ! A row for a certain orbital of a certain atom is as long, as the 433 ! nr. of columns determined previously for the atom. 434 rowpnt(1) = 1 435 ind = 1 436 do iAt1 = 1, nAtom 437 nOrb1 = orb%nOrbAtom(iAt1) 438 do iOrb1 = 1, nOrb1 439 rowpnt(ind+iOrb1) = rowpnt(ind) + iOrb1 * nColAtom(iAt1) 440 end do 441 ind = ind + nOrb1 442 end do 443 deallocate(nColAtom) 444 445 call create(csr, nrow, nrow, rowpnt(nrow + 1) - 1) 446 csr%rowpnt = rowpnt 447 deallocate(rowpnt) 448 449 ! Nr. of CSR columns already filled 450 allocate(nCols(csr%nRow)) 451 nCols(:) = 0 452 ! Index of the nonzero blocks 453 allocate(iNonZero(nAtom)) 454 455 ! Loop over all atoms (over all atomic columns in the rectangular picture) 456 lpAt1: do iAt1 = 1, nAtom 457 iCol = iAtomStart(iAt1) 458 nOrb1 = orb%nOrbAtom(iAt1) ! Width of the atomic column 459 460 ! Mark nonzero blocks in the folded matrix for the current atom 461 nNonZero = 0 462 zero(:) = .true. 463 do iNeigh = 0, nNeighbor(iAt1) 464 iAt2f = img2CentCell(iNeighbor(iNeigh,iAt1)) 465 if (zero(iAt2f)) then 466 zero(iAt2f) = .false. 467 nNonZero = nNonZero + 1 468 iNonzero(nNonzero) = iAt2f 469 end if 470 end do 471 472 ! Initialise CSR internal pointers according the non-zero blocks 473 lpNonZero: do ii = 1, nNonZero 474 iAt2f = iNonZero(ii) 475 nOrb2 = orb%nOrbAtom(iAt2f) 476 iRow = iAtomStart(iAt2f) 477 478 ! Correspond rows of the current atomic block column to the appropriate 479 ! partial rows in the lower triangle of the CSR matrix. 480 do iOrb2 = 0, nOrb2 - 1 481 jj = iRow + iOrb2 482 ind = csr%rowpnt(jj) + nCols(jj) 483 do iOrb1 = 0, nOrb1 - 1 484 csr%colind(ind+iOrb1) = iCol + iOrb1 485 end do 486 nCols(jj) = nCols(jj) + nOrb1 487 end do 488 489 ! Correspond the columns of the current block to appropriate 490 ! partial rows in the upper triangle of the CSR matrix. 491 if (iAt2f /= iAt1) then 492 do iOrb1 = 0, nOrb1 - 1 493 jj = iCol + iOrb1 494 ind = csr%rowpnt(jj) + nCols(jj) 495 do iOrb2 = 0, nOrb2 - 1 496 csr%colind(ind+iOrb2) = iRow + iOrb2 497 end do 498 nCols(jj) = nCols(jj) + nOrb2 499 end do 500 end if 501 end do lpNonZero 502 end do lpAt1 503 504 deallocate(zero) 505 deallocate(nCols) 506 deallocate(iNonZero) 507 508 end subroutine createEquivCSR_cplx 509 510 511 512 !> Folds the internal sparse formatted matrix to the compressed sparse row format (complex 513 !> version). 514 !> 515 !> Note: In the resulting CSR format both triangles of the matrix are set. 516 !> Caveat The routine should only applied on CSR matrixes created by the createEquiv_cplx 517 !> subroutine, since it exploits the ordering of the column indexes. 518 subroutine foldToCSR_cplx(csr, sparse, kPoint, iAtomStart, iPair, iNeighbor, nNeighbor,& 519 & img2CentCell, iCellVec, cellVec, orb) 520 521 !> Resulting CSR matrix. 522 type(z_CSR), intent(inout) :: csr 523 524 !> The internal sparse matrix to fold. 525 real(dp), intent(in) :: sparse(:) 526 527 !> Current k-point 528 real(dp), intent(in) :: kPoint(:) 529 530 !> Starting positions of the atoms in the square matrix 531 integer, intent(in) :: iAtomStart(:) 532 533 !> iPair Starting position of atom-neighbor interaction in the sparse matrix. 534 integer, intent(in) :: iPair(0:,:) 535 536 !> iNeighbor Index of neighbors. 537 integer, intent(in) :: iNeighbor(0:,:) 538 539 !> nNeighbor Number of neighbors. 540 integer, intent(in) :: nNeighbor(:) 541 542 !> img2CentCell Image of the atoms in the central cell. 543 integer, intent(in) :: img2CentCell(:) 544 545 !> Index for cells atomic images are sited in 546 integer, intent(in) :: iCellVec(:) 547 548 !> vectors to crystal unit cells 549 real(dp), intent(in) :: cellVec(:,:) 550 551 !> orb Orbitals in the system. 552 type(TOrbitals), intent(in) :: orb 553 554 integer, allocatable :: nCols(:) 555 complex(dp), allocatable :: tmpCol(:,:), phases(:) 556 integer :: nAtom, nCellVec 557 integer :: iOrb1, nOrb1, nOrb2, iAt1, iAt2, iAt2f, iNeigh 558 integer :: iRow, iCol, iVal, ind 559 integer :: ii, jj, kk 560 561 nAtom = size(orb%nOrbAtom) 562 nCellVec = size(cellVec, dim=2) 563 564 ! Create necessary phase factors 565 allocate(phases(nCellVec)) 566 do ii = 1, nCellVec 567 phases(ii) = exp(2.0_dp * pi * (0.0_dp, 1.0_dp) * dot_product(kPoint, cellVec(:,ii))) 568 end do 569 570 allocate(tmpCol(csr%nRow, orb%mOrb)) ! One atomic block column. 571 lpAt1: do iAt1 = 1, nAtom 572 ! Unfold the columns for the current atom from the sparse matrix. 573 nOrb1 = orb%nOrbAtom(iAt1) 574 tmpCol(:,:) = 0.0_dp 575 do iNeigh = 0, nNeighbor(iAt1) 576 ind = iPair(iNeigh,iAt1) + 1 577 iAt2 = iNeighbor(iNeigh,iAt1) 578 iAt2f = img2CentCell(iAt2) 579 nOrb2 = orb%nOrbAtom(iAt2f) 580 iRow = iAtomStart(iAt2f) 581 tmpCol(iRow:iRow+nOrb2-1, 1:nOrb1) = tmpCol(iRow:iRow+nOrb2-1, 1:nOrb1)& 582 & + phases(iCellVec(iAt2)) * reshape(sparse(ind:ind+nOrb2*nOrb1-1), (/ nOrb2, nOrb1 /)) 583 end do 584 585 ! Copy every column into the appropriate row in the upper triangle of 586 ! the CSR matrix (the copy is masked by the sparsity structure stored 587 ! in the CSR matrix) 588 nOrb1 = orb%nOrbAtom(iAt1) 589 iCol = iAtomStart(iAt1) 590 do iOrb1 = 1, nOrb1 591 ii = csr%rowpnt(iCol + iOrb1 - 1) 592 jj = csr%rowpnt(iCol + iOrb1) - 1 593 do kk=ii,jj 594 csr%nzval(kk) = conjg(tmpCol(csr%colind(kk), iOrb1)) 595 enddo 596 end do 597 end do lpAt1 598 deallocate(tmpCol) 599 deallocate(phases) 600 601 ! Fill the lower triangle of the CSR matrix 602 allocate(nCols(csr%nRow)) 603 nCols(:) = 0 604 do iRow = 1, csr%nRow 605 ! Starting from the tail of the matrix row 606 iVal = csr%rowpnt(iRow + 1) - 1 607 iCol = csr%colind(iVal) 608 do while (iVal >= csr%rowpnt(iRow) .and. iCol > iRow) 609 csr%nzval(csr%rowpnt(iCol)+nCols(iCol)) = conjg(csr%nzval(iVal)) 610 nCols(iCol) = nCols(iCol) + 1 611 iVal = iVal - 1 612 iCol = csr%colind(iVal) 613 end do 614 end do 615 deallocate(nCols) 616 617 end subroutine foldToCSR_cplx 618 619 620 621 !> Unfolds a matrix from the CSR form into the internal sparse representation (real version). 622 !> 623 !> Note: The CSR matrix must be hermitian. The unfolded matrix is added to the passed sparse 624 !> matrix. 625 subroutine unfoldFromCSR_cplx(sparse, csr, kPoint, kWeight, iAtomStart, iPair, iNeighbor,& 626 & nNeighbor, img2CentCell, iCellVec, cellVec, orb) 627 628 !> sparse Updated sparse matrix. 629 real(dp), intent(inout) :: sparse(:) 630 631 !> csr CSR matrix 632 type(z_CSR), intent(in) :: csr 633 634 !> kPoint K-point for the unfolding. 635 real(dp), intent(in) :: kPoint(:) 636 637 !> kWeight Weight of the K-point for the unfolding. 638 real(dp), intent(in) :: kWeight 639 640 !> iAtomStart Starting positions of the atoms in the square matrix 641 integer, intent(in) :: iAtomStart(:) 642 643 !> iPair Starting position of atom-neighbor interaction in the sparse matrix. 644 integer, intent(in) :: iPair(0:,:) 645 646 !> iNeighbor Index of neighbors 647 integer, intent(in) :: iNeighbor(0:,:) 648 649 !> nNeighbor Number of neighbors 650 integer, intent(in) :: nNeighbor(:) 651 652 !> img2CentCell Image of the atoms in the central cell. 653 integer, intent(in) :: img2CentCell(:) 654 655 !> iCellVec Index of the cell shifting vector for each atom. 656 integer, intent(in) :: iCellVec(:) 657 658 !> cellVec Cell shifting vectors (relative coordinates) 659 real(dp), intent(in) :: cellVec(:,:) 660 661 !> orb Orbitals in the system. 662 type(TOrbitals), intent(in) :: orb 663 664 complex(dp), allocatable :: tmpCol(:,:), phases(:) 665 integer :: nAtom, nCellVec 666 integer :: iOrb1, nOrb1, nOrb2, iAt1, iAt2, iAt2f, iNeigh 667 integer :: iRow, iCol, ind 668 integer :: ii 669 670 nAtom = size(orb%nOrbAtom) 671 672 @:ASSERT(csr%nRow == iAtomStart(nAtom+1) - 1) 673 @:ASSERT(size(kPoint) == 3) 674 675 allocate(tmpCol(csr%nRow, orb%mOrb)) 676 677 ! Create phase factors 678 nCellVec = size(cellVec, dim=2) 679 allocate(phases(nCellVec)) 680 do ii = 1, nCellVec 681 phases(ii) = exp(-2.0_dp * pi * (0.0_dp, 1.0_dp) * dot_product(kPoint, cellVec(:,ii))) 682 end do 683 684 do iAt1 = 1, nAtom 685 ! Put the rows belonging to a certain atom into the appropriate column 686 ! of the block column. (Matrix is assumed to be hermitian!) 687 tmpCol(:,:) = (0.0_dp, 0.0_dp) 688 nOrb1 = orb%nOrbAtom(iAt1) 689 iCol = iAtomStart(iAt1) 690 do iOrb1 = 1, nOrb1 691 ii = iCol + iOrb1 - 1 692 do ind = csr%rowpnt(ii), csr%rowpnt(ii+1) - 1 693 tmpCol(csr%colind(ind), iOrb1) = conjg(csr%nzval(ind)) 694 ! NOTE: why conjg ?? 695 end do 696 end do 697 698 ! Unfold the block column into the internal sparse format 699 do iNeigh = 0, nNeighbor(iAt1) 700 ind = iPair(iNeigh,iAt1) + 1 701 iAt2 = iNeighbor(iNeigh, iAt1) 702 iAt2f = img2CentCell(iAt2) 703 nOrb2 = orb%nOrbAtom(iAt2f) 704 iRow = iAtomStart(iAt2f) 705 sparse(ind:ind+nOrb2*nOrb1-1) = sparse(ind:ind+nOrb2*nOrb1-1)& 706 & + kWeight * real( phases(iCellVec(iAt2))& 707 & * reshape(tmpCol(iRow:iRow+nOrb2-1, 1:nOrb1),(/nOrb2*nOrb1/) ), dp) 708 709 end do 710 end do 711 712 deallocate(tmpCol) 713 deallocate(phases) 714 715 end subroutine unfoldFromCSR_cplx 716 717 718 719 !> Creates a CSR matrix by cloning the sparsity structure of another CSR matrix (real version). 720 !> 721 !> Note: The elements of the original matrix are not copied. 722 subroutine cloneSparsityMap_real_real(csrOut, csrIn) 723 724 !> New CSR matrix on exit. 725 type(r_CSR), intent(inout) :: csrOut 726 727 !> CSR matrix with the sparsity map to be cloned. 728 type(r_CSR), intent(in) :: csrIn 729 730 call create(csrOut, csrIn%nRow, csrIn%nCol, csrIn%nnz) 731 csrOut%colind = csrIn%colind 732 csrOut%rowpnt = csrIn%rowpnt 733 csrOut%sorted = csrIn%sorted 734 735 end subroutine cloneSparsityMap_real_real 736 737 738 739 !> Creates a CSR matrix by cloning the sparsity structure of another CSR matrix (complex version). 740 !> 741 !> Note: The elements of the original matrix are not copied. 742 subroutine cloneSparsityMap_cplx_cplx(csrOut, csrIn) 743 744 !> New CSR matrix on exit. 745 type(z_CSR), intent(inout) :: csrOut 746 747 !> CSR matrix with the sparsity map to be cloned. 748 type(z_CSR), intent(in) :: csrIn 749 750 call create(csrOut, csrIn%nRow, csrIn%nCol, csrIn%nnz) 751 csrOut%colind = csrIn%colind 752 csrOut%rowpnt = csrIn%rowpnt 753 csrOut%sorted = csrIn%sorted 754 755 end subroutine cloneSparsityMap_cplx_cplx 756 757 758 759 !> Creates a CSR matrix by cloning the sparsity structure of an other CSR matrix (complex 760 !> version). 761 !> 762 !> Note: The elements of the original matrix are not copied. 763 subroutine cloneSparsityMap_real_cplx(csrOut, csrIn) 764 765 !> New CSR matrix on exit. 766 type(r_CSR), intent(inout) :: csrOut 767 768 !> CSR matrix with the sparsity map to be cloned. 769 type(z_CSR), intent(in) :: csrIn 770 771 call create(csrOut, csrIn%nRow, csrIn%nCol, csrIn%nnz) 772 csrOut%colind = csrIn%colind 773 csrOut%rowpnt = csrIn%rowpnt 774 csrOut%sorted = csrIn%sorted 775 776 end subroutine cloneSparsityMap_real_cplx 777 778 779 780 !> Creates a CSR matrix by cloning the sparsity structure of another CSR matrix (complex version). 781 !> 782 !> Note: The elements of the original matrix are not actually copied (only the sparsity pattern). 783 subroutine cloneSparsityMap_cplx_real(csrOut, csrIn) 784 785 !> New CSR matrix on exit. 786 type(z_CSR), intent(inout) :: csrOut 787 788 !> CSR matrix with the sparsity map to be cloned. 789 type(r_CSR), intent(in) :: csrIn 790 791 call create(csrOut, csrIn%nRow, csrIn%nCol, csrIn%nnz) 792 csrOut%colind = csrIn%colind 793 csrOut%rowpnt = csrIn%rowpnt 794 csrOut%sorted = csrIn%sorted 795 796 end subroutine cloneSparsityMap_cplx_real 797 798 799 800 !> Creates an empty CSR matrix containing zero elements (real version). 801 subroutine createEmptyCSR_real(csr) 802 803 !> Empty CSR matrix on exit. 804 type(r_CSR), intent(out) :: csr 805 806 call create(csr, 0, 0, 0) 807 csr%nnz = 0 808 csr%nRow = 0 809 csr%nCol = 0 810 csr%sorted = .false. 811 812 end subroutine createEmptyCSR_real 813 814 815 816 !> Creates an empty CSR matrix containing zero elements (complex version). 817 subroutine createEmptyCSR_cplx(csr) 818 819 !> Empty CSR matrix on exit. 820 type(z_CSR), intent(out) :: csr 821 822 call create(csr, 0, 0, 0) 823 csr%nnz = 0 824 csr%nRow = 0 825 csr%nCol = 0 826 csr%sorted = .false. 827 828 end subroutine createEmptyCSR_cplx 829 830 !> Destroy a real CSR matrix 831 subroutine rdestroy_CSR(mat) 832 833 !> matrix 834 type(r_CSR) :: mat 835 836 if (allocated(mat%nzval)) then 837 call destroy(mat) 838 end if 839 840 end subroutine rdestroy_CSR 841 842 !> Destroy a real DNS matrix 843 subroutine rdestroy_DNS(mat) 844 845 !> matrix 846 type(r_DNS) :: mat 847 848 if (allocated(mat%val)) then 849 call destroy(mat) 850 end if 851 852 end subroutine rdestroy_DNS 853 854 !> Destroy a complex CSR matrix 855 subroutine zdestroy_CSR(mat) 856 857 !> matrix 858 type(z_CSR) :: mat 859 860 if (allocated(mat%nzval)) then 861 call destroy(mat) 862 end if 863 864 end subroutine zdestroy_CSR 865 866 !> Destroy a complex DNS matrix 867 subroutine zdestroy_DNS(mat) 868 869 !> matrix 870 type(z_DNS) :: mat 871 872 if (allocated(mat%val)) then 873 call destroy(mat) 874 end if 875 876 end subroutine zdestroy_DNS 877 878 879end module dftbp_matconv 880