1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8! This code segment has been fully created by: 9! Nick Papior Andersen, 2013, nickpapior@gmail.com 10! Please conctact the author, prior to re-using this code. 11 12! ************************************************ 13! * Routines for handling the sparsity pattern. * 14! * We supply routines for initialization and * 15! * broadcasting values. * 16! ************************************************ 17 18 19! In general many of the loops below can be followed by examining this loop: 20 21! ! This loop is across the local rows... 22! do lio = 1 , lnr 23! 24! ! Quickly go past the empty regions... (we have nothing to update) 25! if ( l_ncol(lio) == 0 ) cycle 26! 27! ! obtain the global index of the local orbital. 28! io = index_local_to_global(dit,lio,Node) 29! 30! ! Quickly go past the empty regions... (we have nothing to update) 31! if ( up_ncol(io) == 0 ) cycle 32! 33! ! Do a loop in the local sparsity pattern... 34! ! The local sparsity pattern is more "spread", hence 35! ! we do fewer operations by having this as an outer loop 36! do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 37! 38! jo = UCORB(l_col(lind),nr) 39! 40! ! Now search the update region 41! ! This one must *per definition* have less elements. 42! ! Hence, we can exploit this, and find equivalent 43! ! super-cell orbitals. 44! rind = up_ptr(io) 45! ind = rind + SFIND(up_col(rind+1:rind+up_ncol(io)),jo) 46! if ( ind <= rind ) cycle ! The element does not exist 47! 48! ! Obtain the phase for the orbital ij 49! kx = k(1) * xij(1,pnt(lind)) + & 50! k(2) * xij(2,pnt(lind)) + & 51! k(3) * xij(3,pnt(lind)) 52! 53! ph = fact * cdexp(dcmplx(0._dp,-kx)) 54! 55! ! The integration is: 56! ! \rho = e^{-i.k.R} [ \int (Gf^R-Gf^A)/2 dE + \int Gf^R\Gamma Gf^A dE ] 57! 58! if ( non_Eq ) then 59! 60! ! The integration is: 61! ! \rho = e^{-i.k.R} \int Gf^R\Gamma Gf^A dE 62! kx = aimag( ph*zDu(ind) ) 63! dD(lind) = dD(lind) + kx 64! dE(lind) = dE(lind) + aimag( ph*zEu(ind) ) 65! 66! else 67! 68! ! The fact that we have a SYMMETRIC 69! ! update region makes this *tricky* part easy... 70! rin = up_ptr(jo) 71! rind = rin + SFIND(up_col(rin+1:rin+up_ncol(jo)),io) 72! ! We do a check, just to be sure... 73! if ( rind <= rin ) & 74! call die('ERROR: Conjugated symmetrization point does not exist') 75! 76! ! The integration is: 77! ! \rho = e^{-i.k.R} \int (Gf^R-Gf^A)/2 dE 78! 79! dD(lind) = dD(lind) + aimag( ph*(zDu(ind) - conjg(zDu(rind))) ) 80! dE(lind) = dE(lind) + aimag( ph*(zEu(ind) - conjg(zEu(rind))) ) 81! 82! end if 83! 84! end do 85! end do 86 87 88module m_ts_dm_update 89 90 use precision, only : dp 91 use geom_helper, only : UCORB 92 use intrinsic_missing, only : SFIND 93 94 implicit none 95 96 private :: dp 97 private :: ucorb, sfind 98 99contains 100 101 ! *** 102 ! The following scheme should be followed: 103 ! add_*_DM routines are constructed to be able to handle different schemes 104 ! They are called after each k-point and thus the arguments are the following: 105 ! 1. the local update sparsity pattern 106 ! 2. the global update sparsity pattern (suffix 'u') 107 ! The k-point routine is constructed to handle three different methods of doing 108 ! the weighting which is performed here. 109 ! Commonly the routines should accept input matrices 110 ! with the _correct_ sign. 111 ! I.e. on entry non_eq should have + and eq should have - 112 113 subroutine add_k_DM(spDM,spuDM,D_dim2, spEDM, spuEDM, E_dim2, & 114 n_s,sc_off,k, non_Eq) 115 116 use class_OrbitalDistribution 117 use class_Sparsity 118 use class_dSpData2D 119 use class_zSpData2D 120 121! ********************* 122! * INPUT variables * 123! ********************* 124 ! The local integrated sparsity arrays 125 type(dSpData2D), intent(inout) :: spDM, spEDM 126 ! The current k-point global sparsity arrays 127 type(zSpData2D), intent(inout) :: spuDM, spuEDM 128 ! current update region of last dimension 129 integer, intent(in) :: D_dim2, E_dim2 130 ! The k-point 131 real(dp), intent(in) :: k(3) 132 ! The supercell offsets 133 integer, intent(in) :: n_s 134 real(dp), intent(in) :: sc_off(3,0:n_s-1) 135 logical, intent(in) :: non_Eq 136 137 ! Arrays needed for looping the sparsity 138 type(OrbitalDistribution), pointer :: dit 139 type(Sparsity), pointer :: l_s, up_s 140 integer, pointer :: l_ncol(:) , l_ptr(:) , l_col(:) 141 integer, pointer :: up_ncol(:), up_ptr(:), up_col(:) 142 real(dp), pointer :: dD(:,:) , dE(:,:) 143 complex(dp), pointer :: zDu(:,:), zEu(:,:) 144 integer :: lnr, lio, lind, io, ind, nr, jo 145 integer :: rin, rind 146 logical :: hasEDM 147 complex(dp) :: ph(0:n_s-1) 148 149 if ( (.not. initialized(spDM)) .or. (.not. initialized(spuDM)) ) return 150 151 hasEDM = initialized(spEDM) 152 153 ! get the distribution 154 dit => dist(spDM) 155 156 l_s => spar(spDM) 157 call attach(l_s ,n_col=l_ncol,list_ptr=l_ptr,list_col=l_col, & 158 nrows=lnr,nrows_g=nr) 159 dD => val(spDM) 160 if ( hasEDM ) dE => val(spEDM) 161 162 up_s => spar(spuDM) 163 call attach(up_s,n_col=up_ncol,list_ptr=up_ptr,list_col=up_col) 164 zDu => val(spuDM) 165 if ( hasEDM ) zEu => val(spuEDM) 166 167 if ( size(zDu,2) < D_dim2 .or. size(dD,2) < D_dim2 ) then 168 call die('add_k_DM: Error in code') 169 end if 170 171 if ( hasEDM ) then 172 if ( size(zEu,2) < E_dim2 .or. size(dE,2) < E_dim2 ) then 173 call die('add_k_DM: Error in code') 174 end if 175 end if 176 177 ! Remember that this is a sparsity pattern which contains 178 ! a subset of the SIESTA pattern. 179 180 if ( nr /= nrows(up_s) ) call die('The sparsity format is not as & 181 &expected k-DM.') 182 183 do jo = 0 , n_s - 1 184 ph(jo) = cdexp(dcmplx(0._dp, -& 185 k(1) * sc_off(1,jo) - & 186 k(2) * sc_off(2,jo) - & 187 k(3) * sc_off(3,jo))) 188 end do 189 190 if ( non_Eq ) then 191 192 if ( hasEDM ) then 193 194! No data race will occur, sparsity pattern only tranversed once 195!$OMP parallel do default(shared), & 196!$OMP&private(lio,io,lind,jo,rind,ind) 197 do lio = 1 , lnr 198 199 if ( l_ncol(lio) /= 0 ) then 200 io = index_local_to_global(dit,lio) 201 if ( up_ncol(io) /= 0 ) then 202 203 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 204 205 jo = UCORB(l_col(lind),nr) 206 207 rind = up_ptr(io) 208 ind = rind + SFIND(up_col(rind+1:rind+up_ncol(io)),jo) 209 if ( ind <= rind ) cycle ! The element does not exist 210 211 jo = (l_col(lind)-1) / nr 212 213 ! The integration is: 214 ! \rho = e^{-i.k.R} \int Gf^R\Gamma Gf^A dE 215 dD(lind,1:D_dim2) = dD(lind,1:D_dim2) + & 216 real( ph(jo)*zDu(ind,1:D_dim2) ,dp) 217 dE(lind,1:E_dim2) = dE(lind,1:E_dim2) + & 218 real( ph(jo)*zEu(ind,1:E_dim2) ,dp) 219 220 end do 221 222 end if 223 end if 224 225 end do 226!$OMP end parallel do 227 228 else 229 230!$OMP parallel do default(shared), & 231!$OMP&private(lio,io,lind,jo,rind,ind) 232 do lio = 1 , lnr 233 234 if ( l_ncol(lio) /= 0 ) then 235 io = index_local_to_global(dit,lio) 236 if ( up_ncol(io) /= 0 ) then 237 238 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 239 240 jo = UCORB(l_col(lind),nr) 241 242 rind = up_ptr(io) 243 ind = rind + SFIND(up_col(rind+1:rind+up_ncol(io)),jo) 244 if ( ind <= rind ) cycle ! The element does not exist 245 246 jo = (l_col(lind)-1) / nr 247 248 dD(lind,1:D_dim2) = dD(lind,1:D_dim2) + & 249 real( ph(jo)*zDu(ind,1:D_dim2) ,dp) 250 251 end do 252 253 end if 254 end if 255 256 end do 257!$OMP end parallel do 258 259 end if 260 261 else ! non_eq == .false. 262 263 if ( hasEDM ) then 264 265!$OMP parallel do default(shared), & 266!$OMP&private(lio,io,lind,jo,rin,rind,ind) 267 do lio = 1 , lnr 268 269 if ( l_ncol(lio) /= 0 ) then 270 io = index_local_to_global(dit,lio) 271 if ( up_ncol(io) /= 0 ) then 272 273 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 274 275 jo = UCORB(l_col(lind),nr) 276 277 rind = up_ptr(io) 278 ind = rind + SFIND(up_col(rind+1:rind+up_ncol(io)),jo) 279 if ( ind <= rind ) cycle ! The element does not exist 280 281 rin = up_ptr(jo) 282 rind = rin + SFIND(up_col(rin+1:rin+up_ncol(jo)),io) 283#ifndef TS_NOCHECKS 284 if ( rind <= rin ) & 285 call die('ERROR: Conjugated symmetrization point does not exist') 286#endif 287 288 jo = (l_col(lind)-1) / nr 289 290 ! The integration is (- from outside) 291 ! \rho = -\Im e^{-i.k.R} \int (Gf^R-Gf^A) dE 292 dD(lind,1:D_dim2) = dD(lind,1:D_dim2) + & 293 aimag( ph(jo)*(zDu(ind,1:D_dim2) - conjg(zDu(rind,1:D_dim2))) ) 294 dE(lind,1:E_dim2) = dE(lind,1:E_dim2) + & 295 aimag( ph(jo)*(zEu(ind,1:E_dim2) - conjg(zEu(rind,1:E_dim2))) ) 296 297 end do 298 299 end if 300 end if 301 302 end do 303!$OMP end parallel do 304 305 else 306 307!$OMP parallel do default(shared), & 308!$OMP&private(lio,io,lind,jo,rin,rind,ind) 309 do lio = 1 , lnr 310 311 if ( l_ncol(lio) /= 0 ) then 312 io = index_local_to_global(dit,lio) 313 if ( up_ncol(io) /= 0 ) then 314 315 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 316 317 jo = UCORB(l_col(lind),nr) 318 319 rind = up_ptr(io) 320 ind = rind + SFIND(up_col(rind+1:rind+up_ncol(io)),jo) 321 if ( ind <= rind ) cycle ! The element does not exist 322 323 rin = up_ptr(jo) 324 rind = rin + SFIND(up_col(rin+1:rin+up_ncol(jo)),io) 325#ifndef TS_NOCHECKS 326 if ( rind <= rin ) & 327 call die('ERROR: Conjugated symmetrization point does not exist') 328#endif 329 330 jo = (l_col(lind)-1) / nr 331 332 ! - from outside 333 dD(lind,1:D_dim2) = dD(lind,1:D_dim2) + & 334 aimag( ph(jo)*(zDu(ind,1:D_dim2) - conjg(zDu(rind,1:D_dim2))) ) 335 336 end do 337 338 end if 339 end if 340 341 end do 342!$OMP end parallel do 343 344 end if 345 346 end if 347 348 end subroutine add_k_DM 349 350 subroutine add_Gamma_DM(spDM,spuDM,D_dim2,spEDM,spuEDM,E_dim2) 351 352 use class_OrbitalDistribution 353 use class_Sparsity 354 use class_dSpData2D 355 356! ********************* 357! * INPUT variables * 358! ********************* 359 ! The local integrated sparsity arrays 360 type(dSpData2D), intent(inout) :: spDM 361 ! The current Gamma-point global sparsity arrays 362 type(dSpData2D), intent(inout) :: spuDM 363 integer, intent(in) :: D_dim2 364 ! The local integrated sparsity arrays 365 type(dSpData2D), intent(inout) :: spEDM 366 ! The current Gamma-point global sparsity arrays 367 type(dSpData2D), intent(inout) :: spuEDM 368 integer, intent(in) :: E_dim2 369 370 ! Arrays needed for looping the sparsity 371 type(OrbitalDistribution), pointer :: dit 372 type(Sparsity), pointer :: l_s, up_s 373 integer, pointer :: l_ncol(:) , l_ptr(:) , l_col(:) 374 integer, pointer :: up_ncol(:), up_ptr(:), up_col(:) 375 real(dp), pointer :: dD(:,:), dE(:,:) 376 real(dp), pointer :: dDu(:,:), dEu(:,:) 377 integer :: lnr, lio, lind, io, ind, nr, jo, rind 378 logical :: hasEDM 379 380 if ( (.not. initialized(spDM)) .or. (.not. initialized(spuDM)) ) return 381 382 hasEDM = initialized(spEDM) 383 384 ! get distribution 385 dit => dist(spDM) 386 387 l_s => spar(spDM) 388 call attach(l_s ,n_col=l_ncol,list_ptr=l_ptr,list_col=l_col, & 389 nrows=lnr,nrows_g=nr) 390 dD => val(spDM) 391 if ( hasEDM ) dE => val(spEDM) 392 393 up_s => spar(spuDM) 394 call attach(up_s,n_col=up_ncol,list_ptr=up_ptr,list_col=up_col) 395 dDu => val(spuDM) 396 if ( hasEDM ) dEu => val(spuEDM) 397 398 if ( size(dDu,2) < D_dim2 .or. size(dD,2) < D_dim2 ) then 399 call die('add_Gamma_DM: Error in code') 400 end if 401 402 if ( hasEDM ) then 403 if ( size(dEu,2) < E_dim2 .or. size(dE,2) < E_dim2 ) then 404 call die('add_Gamma_DM: Error in code') 405 end if 406 end if 407 408 ! Remember that this is a sparsity pattern which contains 409 ! a subset of the SIESTA pattern. 410 411 if ( nr /= nrows(up_s) ) call die('The sparsity format is not as & 412 &expected G-DM.') 413 414 if ( hasEDM ) then 415 416!$OMP parallel do default(shared), & 417!$OMP&private(lio,io,lind,jo,rind,ind) 418 do lio = 1 , lnr 419 420 if ( l_ncol(lio) /= 0 ) then 421 io = index_local_to_global(dit,lio) 422 if ( up_ncol(io) /= 0 ) then 423 424 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 425 426 ! we might still have SIESTA-non-Gamma 427 jo = ucorb(l_col(lind),nr) 428 429 ! This sparsity pattern is in UC 430 rind = up_ptr(io) 431 ind = rind + SFIND(up_col(rind+1:rind+up_ncol(io)),jo) 432 if ( ind <= rind ) cycle ! The element does not exist 433 434 dD(lind,1:D_dim2) = dD(lind,1:D_dim2) + dDu(ind,1:D_dim2) 435 dE(lind,1:E_dim2) = dE(lind,1:E_dim2) + dEu(ind,1:E_dim2) 436 437 end do 438 439 end if 440 end if 441 442 end do 443!$OMP end parallel do 444 445 else 446 447!$OMP parallel do default(shared), & 448!$OMP&private(lio,io,lind,jo,rind,ind) 449 do lio = 1 , lnr 450 451 if ( l_ncol(lio) /= 0 ) then 452 io = index_local_to_global(dit,lio) 453 if ( up_ncol(io) /= 0 ) then 454 455 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 456 457 jo = ucorb(l_col(lind),nr) 458 459 rind = up_ptr(io) 460 ind = rind + SFIND(up_col(rind+1:rind+up_ncol(io)),jo) 461 if ( ind <= rind ) cycle ! The element does not exist 462 463 dD(lind,1:D_dim2) = dD(lind,1:D_dim2) + dDu(ind,1:D_dim2) 464 465 end do 466 467 end if 468 end if 469 470 end do 471!$OMP end parallel do 472 473 end if 474 475 end subroutine add_Gamma_DM 476 477 subroutine update_DM(dit,sp,n_nzs,DM, spDM, Ef, & 478 EDM, spEDM, ipnt, UpSpGlobal) 479 480 use class_OrbitalDistribution 481 use class_Sparsity 482 use class_dSpData2D 483 use class_iSpData1D 484 485! ********************* 486! * INPUT variables * 487! ********************* 488 type(OrbitalDistribution), intent(inout) :: dit 489 type(Sparsity), intent(inout) :: sp 490 ! Size of the sparsity arrays 491 integer, intent(in) :: n_nzs 492 ! Sparse DM-arrays (local) 493 real(dp), intent(inout) :: DM(n_nzs) 494 ! Updated sparsity arrays (they contain the current integration) 495 type(dSpData2D), intent(inout) :: spDM 496 ! fermi-level, we shift the energy density matrix back 497 real(dp), intent(in) :: Ef 498 ! Sparse energy-DM-arrays (local) 499 real(dp), intent(inout) :: EDM(n_nzs) 500 ! Updated sparsity arrays (they contain the current integration) 501 type(dSpData2D), intent(inout) :: spEDM 502 ! I.e. a pointer from the local update sparsity to the local sparsity 503 type(iSpData1D), intent(in), optional :: ipnt 504 ! Whether the update sparsity pattern is a global update sparsity pattern 505 logical, intent(in), optional :: UpSpGlobal 506 507 ! Arrays needed for looping the sparsity 508 type(Sparsity), pointer :: s 509 integer, pointer :: l_ncol(:), l_ptr(:), l_col(:) 510 integer, pointer :: lup_ncol(:), lup_ptr(:), lup_col(:) 511 integer, pointer :: pnt(:) 512 real(dp), pointer :: dD(:,:), dE(:,:) 513 integer :: lnr, nr, uind, lio, io, lind, ind, ljo, jo 514 logical :: hasipnt, hasEDM, lUpSpGlobal 515 516 call attach(sp, n_col=l_ncol,list_ptr=l_ptr,list_col=l_col, & 517 nrows=lnr,nrows_g=nr) 518 s => spar(spDM) 519 call attach(s, n_col=lup_ncol,list_ptr=lup_ptr,list_col=lup_col) 520 dD => val(spDM) 521 522 hasEDM = initialized(spEDM) 523 524 if ( hasEDM ) then 525 526 dE => val(spEDM) 527 528 ! the actual size of the shift 529 jo = nnzs(spDM) 530 531 ! As we have shifted the fermi-level up to 0, we need to shift the 532 ! energy-density matrix back 533 call daxpy(jo,Ef,dD(1,1),1,dE(1,1),1) 534 535 end if 536 537 ! We have that the update sparsity pattern is in local 538 ! form. 539 ! this means that sp == s (besides the non-update objects) 540 ! Hence we don't need to utilize index_local_to_global 541 542 543 hasipnt = present(ipnt) 544 if ( hasipnt ) hasipnt = initialized(ipnt) 545 546 lUpSpGlobal = .false. 547 if ( present(UpSpGlobal) ) lUpSpGlobal = UpSpGlobal 548 549 if ( .not. lUpSpGlobal ) then 550 551 if ( lnr /= nrows(s) ) & 552 call die('The sparsity format is not as expected u-DM.') 553 554 if ( hasipnt ) then 555 556 ! The pointer 557 pnt => val(ipnt) 558 559 ! This loop is across the local rows... 560! No data race will occur 561!$OMP parallel do default(shared), & 562!$OMP&private(io,uind,ind) 563 do io = 1 , lnr 564 565 ! Quickly go past the empty regions... (we have nothing to update) 566 if ( lup_ncol(io) /= 0 ) then 567 568 ! Do a loop in the local update sparsity pattern... 569 ! The local sparsity pattern is more "spread", hence 570 ! we do fewer operations by having this as an outer loop 571 do uind = lup_ptr(io) + 1 , lup_ptr(io) + lup_ncol(io) 572 573 ind = pnt(uind) 574 575 DM(ind) = DM(ind) + dD(uind,1) 576 if ( hasEDM ) EDM(ind) = EDM(ind) + dE(uind,1) 577 578 end do 579 580 end if 581 end do 582!$OMP end parallel do 583 584 else 585 586 ! This loop is across the local rows... 587! no data race will occur 588!$OMP parallel do default(shared), & 589!$OMP&private(io,uind,jo,ind) 590 do io = 1 , lnr 591 592 ! Quickly go past the empty regions... (we have nothing to update) 593 if ( lup_ncol(io) /= 0 ) then 594 595 ! Do a loop in the local update sparsity pattern... 596 ! The local sparsity pattern is more "spread", hence 597 ! we do fewer operations by having this as an outer loop 598 do uind = lup_ptr(io) + 1 , lup_ptr(io) + lup_ncol(io) 599 600 ! We are dealing with a non-UC sparsity pattern 601 jo = lup_col(uind) 602 603 ! Now we loop across the local region 604 ind = l_ptr(io) + & 605 minloc(abs(l_col(l_ptr(io)+1:l_ptr(io)+l_ncol(io))-jo),1) 606 if ( l_col(ind) /= jo ) then 607 do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io) 608 if ( l_col(ind) /= jo ) cycle 609 exit ! we have the correct ind-value 610 end do 611 if ( l_col(ind) /= jo ) cycle 612 end if 613 614 ! We need to add in case of special weighting... 615 DM(ind) = DM(ind) + dD(uind,1) 616 if ( hasEDM ) EDM(ind) = EDM(ind) + dE(uind,1) 617 618 end do 619 end if 620 end do 621!$OMP end parallel do 622 end if 623 624 else 625 ! We have a global update sparsity pattern 626 627 ! This is the global sparsity pattern 628 ! i.e. we require to call index_local_to_global 629 ! The global sparsity pattern is not in supercell format 630 631 ! This loop is across the local rows... 632! We will never have a data race here (it is on local sparsity pattern) 633!$OMP parallel do default(shared), & 634!$OMP&private(lio,io,lind,jo,ljo,ind) 635 do lio = 1 , lnr 636 637 ! obtain the global index of the local orbital. 638 io = index_local_to_global(dit,lio) 639 640 ! Quickly go past the empty regions... (we have nothing to update) 641 if ( lup_ncol(io) /= 0 ) then 642 643 ! Retrieve pointer index 644 jo = lup_ptr(io) 645 646 ! Do a loop in the local sparsity pattern... 647 ! The local sparsity pattern is more "spread", hence 648 ! we do fewer operations by having this as an outer loop 649 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 650 651 ! we need to compare with the global update sparsity 652 ljo = UCORB(l_col(lind),nr) 653 654 ! Now we loop across the update region 655 ! This one must *per definition* have less elements. 656 ! Hence, we can exploit this, and find equivalent 657 ! super-cell orbitals. 658 ! Ok, this is Gamma (but to be consistent) 659 ind = jo + SFIND(lup_col(jo+1:jo+lup_ncol(io)),ljo) 660 if ( ind > jo ) then 661 DM(lind) = DM(lind) + dD(ind,1) 662 if ( hasEDM ) EDM(lind) = EDM(lind) + dE(ind,1) 663 end if 664 665 end do 666 667 end if 668 end do 669!$OMP end parallel do 670 671 end if 672 673 end subroutine update_DM 674 675 ! This routine will ONLY be called if .not. IsVolt, 676 ! Hence we don't have any sparsity patterns with local sparsity patterns 677 ! that is dealing with this routine (hence we do need the index_local_to_global) 678 subroutine update_zDM(dit,sp,n_nzs,DM,spDM, Ef, & 679 EDM,spEDM, k, n_s, sc_off) 680 use class_OrbitalDistribution 681 use class_Sparsity 682 use class_zSpData2D 683 684 type(OrbitalDistribution), intent(inout) :: dit 685 type(Sparsity), intent(inout) :: sp 686 ! Size of the sparsity arrays 687 integer, intent(in) :: n_nzs 688 ! Sparse DM-arrays (local) 689 real(dp), intent(inout) :: DM(n_nzs), EDM(n_nzs) 690 ! Updated sparsity arrays (they contain the current integration) 691 type(zSpData2D), intent(inout) :: spDM, spEDM 692 ! The fermi level 693 real(dp), intent(in) :: Ef 694 ! The k-point... 695 real(dp), intent(in) :: k(3) 696 ! The supercell offset 697 integer, intent(in) :: n_s 698 real(dp), intent(in) :: sc_off(3,0:n_s-1) 699 700 ! Arrays needed for looping the sparsity 701 type(Sparsity), pointer :: s 702 integer, pointer :: l_ncol(:), l_ptr(:), l_col(:) 703 integer, pointer :: lup_ncol(:), lup_ptr(:), lup_col(:) 704 complex(dp), pointer :: zD(:,:), zE(:,:) 705 complex(dp) :: ph(0:n_s-1) 706 integer :: lio, io, jo, ind, nr 707 integer :: lnr, lind, rin, rind 708 logical :: hasEDM 709 710 call attach(sp, n_col=l_ncol,list_ptr=l_ptr,list_col=l_col, & 711 nrows=lnr,nrows_g=nr) 712 s => spar(spDM) 713 call attach(s, n_col=lup_ncol,list_ptr=lup_ptr,list_col=lup_col) 714 zD => val(spDM) 715 716 hasEDM = initialized(spEDM) 717 718 if ( hasEDM ) then 719 720 zE => val(spEDM) 721 722 ! the actual size of the shift 723 jo = nnzs(spDM) 724 725 ! As we have shifted the fermi-level up to 0, we need to shift the 726 ! energy-density matrix back 727 ph(0) = dcmplx(Ef,0._dp) 728 call zaxpy(jo,ph(0),zD(1,1),1,zE(1,1),1) 729 730 end if 731 732 ! Remember that this is a sparsity pattern which contains 733 ! a subset of the SIESTA pattern (but still the global sparsity pattern) 734 735 if ( nr /= nrows(s) ) call die('The sparsity format is not as & 736 &expected uz-DM.') 737 738 do jo = 0 , n_s - 1 739 ph(jo) = cdexp(dcmplx(0._dp, -& 740 k(1) * sc_off(1,jo) - & 741 k(2) * sc_off(2,jo) - & 742 k(3) * sc_off(3,jo))) 743 end do 744 745 ! This loop is across the local rows... 746! No data race will occur 747!$OMP parallel do default(shared), & 748!$OMP&private(lio,io,lind,jo,rind,ind,rin) 749 do lio = 1 , lnr 750 751 ! obtain the global index of the local orbital. 752 io = index_local_to_global(dit,lio) 753 754 ! Quickly go past the empty regions... (we have nothing to update) 755 if ( lup_ncol(io) /= 0 ) then 756 757 ! Do a loop in the local sparsity pattern... 758 ! The local sparsity pattern is more "spread", hence 759 ! we do fewer operations by having this as an outer loop 760 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 761 762 jo = UCORB(l_col(lind),nr) 763 764 ! Now search the update region 765 ! This one must *per definition* have less elements. 766 ! Hence, we can exploit this, and find equivalent 767 ! super-cell orbitals. 768 rind = lup_ptr(io) 769 ind = rind + SFIND(lup_col(rind+1:rind+lup_ncol(io)),jo) 770 if ( ind <= rind ) cycle ! The element does not exist 771 772 ! The fact that we have a SYMMETRIC 773 ! update region makes this *tricky* part easy... 774 rin = lup_ptr(jo) 775 ! TODO, this REQUIRES that lup_col(:) is sorted 776 rind = rin + SFIND(lup_col(rin+1:rin+lup_ncol(jo)),io) 777#ifndef TS_NOCHECKS 778 ! We do a check, just to be sure... 779 if ( rind <= rin ) & 780 call die('ERROR: Conjugated symmetrization point does not exist') 781#endif 782 783 ! The integration is: (- is from outside) 784 ! \rho = - \Im e^{-i.k.R} [ \int (Gf^R-Gf^A) dE + \int Gf^R\Gamma Gf^A dE ] 785 jo = (l_col(lind)-1) / nr 786 787 DM(lind) = DM(lind) + aimag( ph(jo)*(zD(ind,1) - conjg(zD(rind,1))) ) 788 if ( hasEDM ) & 789 EDM(lind) = EDM(lind) + aimag( ph(jo)*(zE(ind,1) - conjg(zE(rind,1))) ) 790 791 end do 792 793 end if 794 end do 795!$OMP end parallel do 796 797 end subroutine update_zDM 798 799 800 subroutine init_DM(dit,sp,n_nzs,DM,EDM, up_sp, Calc_Forces) 801 ! The DM and EDM equivalent matrices 802 use class_OrbitalDistribution 803 use class_Sparsity 804 805 type(OrbitalDistribution), intent(inout) :: dit 806 type(Sparsity), intent(inout) :: sp 807 ! Size of the sparsity arrays 808 integer, intent(in) :: n_nzs 809 ! Sparse DM-arrays (local) 810 real(dp), intent(inout) :: DM(n_nzs), EDM(n_nzs) 811 ! The updated sparsity arrays... 812 type(Sparsity), intent(inout) :: up_sp 813 logical, intent(in) :: Calc_Forces 814 815 integer, pointer :: l_ncol(:), l_ptr(:), l_col(:) 816 integer, pointer :: lup_ncol(:), lup_ptr(:), lup_col(:) 817 integer :: lnr, lio, lind, io, jo, ind, nr, col 818 819 call attach(sp, n_col=l_ncol,list_ptr=l_ptr,list_col=l_col, & 820 nrows=lnr,nrows_g=nr) 821 call attach(up_sp, n_col=lup_ncol,list_ptr=lup_ptr,list_col=lup_col) 822 823 ! Remember that this is a sparsity pattern which contains 824 ! a subset of the SIESTA pattern. 825 if ( nr /= nrows(up_sp) ) call die('The sparsity format is not as & 826 &expected i-DM.') 827 828 if ( Calc_Forces ) then 829 830 ! This loop is across the local rows... 831!$OMP parallel do default(shared), & 832!$OMP&private(lio,io,jo,col,ind,lind) 833 do lio = 1 , lnr 834 835 ! obtain the global index of the local orbital. 836 io = index_local_to_global(dit,lio) 837 838 ! Quickly go past the empty regions... (we have nothing to update) 839 if ( lup_ncol(io) /= 0 ) then 840 841 jo = lup_ptr(io) 842 843 ! Do a loop in the local sparsity pattern... 844 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 845 846 ! Search for unit-cell entry in update sparsity 847 ! pattern (UC) 848 col = UCORB(l_col(lind),nr) 849 ind = SFIND(lup_col(jo+1:jo+lup_ncol(io)),col) 850 if ( ind > 0 ) then 851 DM (lind) = 0._dp 852 EDM(lind) = 0._dp 853 end if 854 855 end do 856 end if 857 end do 858!$OMP end parallel do 859 860 else 861!$OMP parallel do default(shared), & 862!$OMP&private(lio,io,jo,col,ind,lind) 863 do lio = 1 , lnr 864 io = index_local_to_global(dit,lio) 865 if ( lup_ncol(io) /= 0 ) then 866 jo = lup_ptr(io) 867 do lind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio) 868 col = UCORB(l_col(lind),nr) 869 ind = SFIND(lup_col(jo+1:jo+lup_ncol(io)),col) 870 if ( ind > 0 ) DM (lind) = 0._dp 871 end do 872 end if 873 end do 874!$OMP end parallel do 875 end if 876 877 end subroutine init_DM 878 879end module m_ts_dm_update 880