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#:include "common.fypp" 8 9!> Interface to LIBNEGF for DFTB+ 10module negf_int 11 use libnegf_vars 12 use libnegf, only : convertcurrent, eovh, getel, lnParams, pass_DM, Tnegf, unit 13#:if WITH_MPI 14 use libnegf, only : negf_mpi_init 15#:endif 16 use libnegf, only : z_CSR, READ_SGF, COMP_SGF, COMPSAVE_SGF 17 use libnegf, only : associate_lead_currents, associate_ldos, associate_transmission 18 use libnegf, only : associate_current, compute_current, compute_density_dft, compute_ldos 19 use libnegf, only : create, create_scratch, destroy, set_readoldDMsgf 20 use libnegf, only : destroy_matrices, destroy_negf, get_params, init_contacts, init_ldos 21 use libnegf, only : init_negf, init_structure, pass_hs, set_bp_dephasing 22 use libnegf, only : set_drop, set_elph_block_dephasing, set_elph_dephasing, set_elph_s_dephasing 23 use libnegf, only : set_ldos_indexes, set_params, set_scratch, writememinfo, writepeakinfo 24 use libnegf, only : printcsr 25 use dftbp_accuracy 26 use dftbp_constants 27 use dftbp_matconv 28 use dftbp_sparse2dense 29 use dftbp_densedescr 30 use dftbp_commontypes, only : TOrbitals 31 use dftbp_formatout 32 use dftbp_globalenv 33 use dftbp_message 34 use dftbp_elecsolvertypes, only : electronicSolverTypes 35 use dftbp_linkedlist 36 use dftbp_periodic 37 use dftbp_assert 38#:if WITH_MPI 39 use dftbp_mpifx 40#:endif 41 implicit none 42 private 43 44 type(TNegf), target, public :: negf 45 ! Workaround: ifort 17, ifort 16 46 ! Passing negf for pointer dummy arguments fails despite target attribute, so pointer is needed 47 type(TNegf), pointer :: pNegf 48 49 !> general library initializations 50 public :: negf_init 51 52 !> passing structure parameters 53 public :: negf_init_str 54 55 public :: negf_init_dephasing 56 public :: negf_init_elph 57 public :: negf_destroy 58 59 !> wrapped functions passing dftb matrices. Needed for parallel 60 public :: calcdensity_green 61 62 !> Calculate the energy weighted density from the GF 63 public :: calcEdensity_green 64 65 !> Calculate the total current 66 public :: calc_current 67 68 !> calculate local currents 69 public :: local_currents 70 71 !> interface csr matrices. The pattering must be predefined using negf_init_csr 72 public :: negf_init_csr 73 74 !> compressed sparse row hamiltonian 75 type(z_CSR), target :: csrHam 76 type(Z_CSR), pointer :: pCsrHam => null() 77 78 !> compressed sparse row overlap 79 type(z_CSR), target :: csrOver 80 type(Z_CSR), pointer :: pCsrOver => null() 81 82 !> non wrapped direct calls 83 private :: negf_density, negf_current, negf_ldos 84 85 contains 86 87 !> Init gDFTB environment and variables 88 !> 89 !> Note: mpicomm should be the global commworld here 90#:if WITH_MPI 91 subroutine negf_init(transpar, greendens, tundos, mpicomm, tempElec, solver) 92 93 !> MPI communicator 94 type(mpifx_comm), intent(in) :: mpicomm 95#:else 96 subroutine negf_init(transpar, greendens, tundos, tempElec, solver) 97#:endif 98 99 !> Parameters for the transport calculation 100 Type(TTranspar), intent(in) :: transpar 101 102 !> Parameters for the Green's function calculation 103 Type(TNEGFGreenDensInfo), intent(in) :: greendens 104 105 !> parameters for tuneling and density of states evaluation 106 Type(TNEGFTunDos), intent(in) :: tundos 107 108 !> Electronic temperature 109 real(dp), intent(in) :: tempElec 110 111 !> Which solver call is used in the main code 112 integer, intent(in) :: solver 113 114 ! local variables 115 real(dp), allocatable :: pot(:), eFermi(:) 116 integer :: i, l, ncont, nc_vec(1), j, nldos 117 integer, allocatable :: sizes(:) 118 type(lnParams) :: params 119 120 ! Workaround: ifort 16 121 ! Pointer must be set within a subroutine. Initialization at declaration fails. 122 pNegf => negf 123#:if WITH_MPI 124 call negf_mpi_init(mpicomm) 125#:endif 126 127 if (transpar%defined) then 128 ncont = transpar%ncont 129 else 130 ncont = 0 131 endif 132 133 ! ------------------------------------------------------------------------------ 134 ! Set defaults and fill up the parameter structure with them 135 call init_negf(negf) 136 call init_contacts(negf, ncont) 137 call set_scratch(negf, ".") 138 139 if (tIoProc .and. greendens%saveSGF ) then 140 call create_scratch(negf) 141 end if 142 143 call get_params(negf, params) 144 145 ! ------------------------------------------------------------------------------ 146 ! This must be different for different initialisations, to be separated 147 ! Higher between transport and greendens is taken, temporary 148 if (tundos%defined .and. greendens%defined) then 149 if (tundos%verbose.gt.greendens%verbose) then 150 params%verbose = tundos%verbose 151 else 152 params%verbose = greendens%verbose 153 endif 154 else 155 if (tundos%defined) then 156 params%verbose = tundos%verbose 157 end if 158 if (greendens%defined) then 159 params%verbose = greendens%verbose 160 end if 161 end if 162 ! ------------------------------------------------------------------------------ 163 ! This parameter is used to set the averall drop threshold in libnegf 164 ! It affects especially transmission that is not accurate more than 165 ! this value. 166 call set_drop(1.0e-20_dp) 167 168 169 ! ------------------------------------------------------------------------------ 170 ! Assign spin degenracy and check consistency between different input blocks 171 if (tundos%defined .and. greendens%defined) then 172 if (tundos%gSpin /= greendens%gSpin) then 173 call error("spin degeneracy is not consistent between different input blocks") 174 else 175 params%g_spin = real(tundos%gSpin) ! Spin degeneracy 176 end if 177 else if (tundos%defined) then 178 params%g_spin = real(tundos%gSpin) ! Spin degeneracy 179 else if (greendens%defined) then 180 params%g_spin = real(greendens%gSpin) ! Spin degeneracy 181 end if 182 183 184 ! Setting contact temperatures 185 do i = 1, ncont 186 if (greendens%defined) then 187 params%kbT_dm(i) = greendens%kbT(i) 188 else 189 params%kbT_dm(i) = tempElec 190 end if 191 if (tundos%defined) then 192 params%kbT_t(i) = tundos%kbT(i) 193 else 194 params%kbT_t(i) = tempElec 195 end if 196 end do 197 198 if (ncont == 0) then 199 if (greendens%defined) then 200 params%kbT_dm(1) = greendens%kbT(1) 201 else 202 params%kbT_dm(1) = tempElec 203 end if 204 end if 205 206 ! Make sure low temperatures (< 10K) converted to 0. 207 ! This avoid numerical problems with contour integration 208 do i = 1, size(params%kbT_dm) 209 if (params%kbT_dm(i) < 3.0e-5_dp) then 210 params%kbT_dm(i) = 0.0_dp 211 end if 212 end do 213 214 ! ------------------------------------------------------------------------------ 215 ! SETTING ELECTROCHEMICAL POTENTIALS INCLUDING BUILT-IN 216 ! ------------------------------------------------------------------------------ 217 ! Fermi level is given by the contacts. If no contacts => no transport, 218 ! Then Fermi is defined by the Green solver 219 if (transpar%defined) then 220 pot = transpar%contacts(1:ncont)%potential 221 eFermi = transpar%contacts(1:ncont)%eFermi(1) 222 do i = 1,ncont 223 ! Built-in potential to equilibrate Fermi levels 224 pot(i) = pot(i) + eFermi(i) - minval(eFermi(1:ncont)) 225 226 ! set parameters for wide band approximations 227 params%FictCont(i) = transpar%contacts(i)%wideBand 228 params%contact_DOS(i) = transpar%contacts(i)%wideBandDOS 229 230 write(stdOut,*) '(negf_init) CONTACT INFO #',i 231 232 if (params%FictCont(i)) then 233 write(stdOut,*) 'FICTITIOUS CONTACT ' 234 write(stdOut,*) 'DOS: ', params%contact_DOS(i) 235 end if 236 write(stdOut,*) 'Temperature (DM): ', params%kbT_dm(i) 237 write(stdOut,*) 'Temperature (Current): ', params%kbT_t(i) 238 write(stdOut,*) 'Potential (with built-in): ', pot(i) 239 write(stdOut,*) 'eFermi: ', eFermi(i) 240 write(stdOut,*) 241 242 enddo 243 244 ! Define electrochemical potentials 245 params%mu(1:ncont) = eFermi(1:ncont) - pot(1:ncont) 246 write(stdOut,*) 'Electro-chemical potentials: ', params%mu(1:ncont) 247 write(stdOut,*) 248 deallocate(pot) 249 250 else 251 params%mu(1) = greendens%oneFermi(1) 252 end if 253 254 ! ------------------------------------------------------------------------------ 255 ! SETTING COUNTOUR INTEGRATION PARAMETERS 256 ! ------------------------------------------------------------------------------ 257 if (greendens%defined) then 258 params%Ec = greendens%enLow ! lowest energy 259 params%Np_n(1:2) = greendens%nP(1:2) ! contour npoints 260 params%n_kt = greendens%nkt ! n*kT for Fermi 261 262 ! Real-axis points. 263 ! Override to 0 if bias is 0.0 264 params%Np_real = 0 265 if (ncont > 0) then 266 if (any(abs(params%mu(2:ncont)-params%mu(1)) > 1.0e-10_dp)) then 267 params%Np_real = greendens%nP(3) ! real axis points 268 end if 269 end if 270 271 ! Setting for Read/Write Surface GFs. 272 ! NOTE: for the moment in tunneling and dos SGF are always 273 ! recomputed because bias may change points and errors are easy 274 275 ! Read G.F. from very first iter 276 if (greendens%readSGF) then 277 params%readOldDM_SGFs = READ_SGF 278 end if 279 ! Compute G.F. at every iteration 280 if (.not.greendens%readSGF .and. .not.greendens%saveSGF) then 281 params%readOldDM_SGFs = COMP_SGF 282 end if 283 ! Default Write on first iter 284 if (.not.greendens%readSGF .and. greendens%saveSGF) then 285 params%readOldDM_SGFs = COMPSAVE_SGF 286 end if 287 288 if(any(params%kbT_dm > 0) .and. greendens%nPoles == 0) then 289 call error("Number of Poles = 0 but T > 0") 290 else 291 params%n_poles = greendens%nPoles 292 end if 293 if(all(params%kbT_dm.eq.0)) then 294 params%n_poles = 0 295 end if 296 297 write(stdOut,*) 'Density Matrix Parameters' 298 if (.not.transpar%defined) then 299 write(stdOut,*) 'Temperature (DM): ', params%kbT_dm(1) 300 write(stdOut,*) 'eFermi: ', params%mu(1) 301 end if 302 write(stdOut,*) 'Contour Points: ', params%Np_n(1:2) 303 write(stdOut,*) 'Number of poles: ', params%N_poles 304 write(stdOut,*) 'Real-axis points: ', params%Np_real(1) 305 if (params%readOldDM_SGFs==0) then 306 write(stdOut,*) 'Read Existing SGFs: Yes ' 307 else 308 write(stdOut,*) 'Read Existing SGFs: No, option ', params%readOldDM_SGFs 309 end if 310 write(stdOut,*) 311 312 end if 313 314 ! ------------------------------------------------------------------------------ 315 ! Setting the delta: priority on Green Solver, if present 316 ! dos_delta is used by libnegf to smoothen T(E) and DOS(E) 317 ! and is currently set in tunneling 318 if (tundos%defined) then 319 params%dos_delta = tundos%broadeningDelta 320 params%delta = tundos%delta ! delta for G.F. 321 end if 322 if (greendens%defined) then 323 params%delta = greendens%delta ! delta for G.F. 324 end if 325 326 ! ------------------------------------------------------------------------------ 327 ! SETTING TRANSMISSION PARAMETERS 328 ! ------------------------------------------------------------------------------ 329 if (tundos%defined) then 330 331 l = size(tundos%ni) 332 params%ni(1:l) = tundos%ni(1:l) 333 params%nf(1:l) = tundos%nf(1:l) 334 335 ! setting of intervals and indices for projected DOS 336 nldos = size(tundos%dosOrbitals) 337 call init_ldos(negf, nldos) 338 do i = 1, nldos 339 call set_ldos_indexes(negf, i, tundos%dosOrbitals(i)%data) 340 end do 341 342 params%Emin = tundos%Emin 343 params%Emax = tundos%Emax 344 params%Estep = tundos%Estep 345 346 ! For the moment tunneling and ldos SGFs are always recomputed 347 params%readOldT_SGFs = COMP_SGF 348 349 endif 350 351 ! Energy conversion only affects output units. 352 ! The library writes energies as (E * negf%eneconv) 353 params%eneconv = Hartree__eV 354 355 if (allocated(sizes)) then 356 deallocate(sizes) 357 end if 358 359 call set_params(negf,params) 360 361 !-------------------------------------------------------------------------- 362 ! DAR begin - negf_init - TransPar to negf 363 !-------------------------------------------------------------------------- 364 if (transpar%defined) then 365 !negf%tNoGeometry = transpar%tNoGeometry 366 negf%tOrthonormal = transpar%tOrthonormal 367 negf%tOrthonormalDevice = transpar%tOrthonormalDevice 368 negf%NumStates = transpar%NumStates 369 negf%tManyBody = transpar%tManyBody 370 negf%tElastic = transpar%tElastic 371 negf%tZeroCurrent = transpar%tZeroCurrent 372 negf%MaxIter = transpar%MaxIter 373 negf%trans%out%tWriteDOS = transpar%tWriteDOS 374 negf%tWrite_ldos = transpar%tWrite_ldos 375 negf%tWrite_negf_params = transpar%tWrite_negf_params 376 negf%trans%out%tDOSwithS = transpar%tDOSwithS 377 negf%cont(:)%name = transpar%contacts(:)%name 378 negf%cont(:)%tWriteSelfEnergy = transpar%contacts(:)%tWriteSelfEnergy 379 negf%cont(:)%tReadSelfEnergy = transpar%contacts(:)%tReadSelfEnergy 380 negf%cont(:)%tWriteSurfaceGF = transpar%contacts(:)%tWriteSurfaceGF 381 negf%cont(:)%tReadSurfaceGF = transpar%contacts(:)%tReadSurfaceGF 382 end if 383 384 ! Defined outside transpar%defined ... HAS TO BE FIXED 385 negf%tDephasingVE = transpar%tDephasingVE 386 negf%tDephasingBP = transpar%tDephasingBP 387 388 389 if((.not.negf%tElastic).and.(.not.negf%tManyBody)) then 390 write(stdOut, *)'Current is not calculated!' 391 call error('Choose "Elastic = Yes" or "ManyBody = Yes"!') 392 end if 393 394 395 end subroutine negf_init 396 397 398 !> Initialise dephasing effects 399 subroutine negf_init_dephasing(tundos) 400 401 !> density of states in tunnel region 402 Type(TNEGFTunDos), intent(in) :: tundos 403 404 if(negf%tDephasingVE) then 405 call negf_init_elph(tundos%elph) 406 end if 407 408 if(negf%tDephasingBP) then 409 call negf_init_bp(tundos%bp) 410 end if 411 412 end subroutine negf_init_dephasing 413 414 415 !> Initialise electron-phonon coupling model 416 subroutine negf_init_elph(elph) 417 418 !> el-ph coupling structure 419 type(TElPh), intent(in) :: elph 420 421 write(stdOut,*) 422 select case(elph%model) 423 case(1) 424 write(stdOut,*) 'Setting local fully diagonal (FD) elastic dephasing model' 425 call set_elph_dephasing(negf, elph%coupling, elph%scba_niter) 426 case(2) 427 write(stdOut,*) 'Setting local block diagonal (BD) elastic dephasing model' 428 call set_elph_block_dephasing(negf, elph%coupling, elph%orbsperatm, elph%scba_niter) 429 case(3) 430 write(stdOut,*) 'Setting overlap mask (OM) block diagonal elastic dephasing model' 431 call set_elph_s_dephasing(negf, elph%coupling, elph%orbsperatm, elph%scba_niter) 432 case default 433 call error("This electron-phonon model is not supported") 434 end select 435 436 end subroutine negf_init_elph 437 438 439 !> Initialise Buttiker Probe dephasing 440 subroutine negf_init_bp(elph) 441 442 !> el-ph coupling structure 443 type(TElPh), intent(in) :: elph 444 445 write(stdOut,*) 446 select case(elph%model) 447 case(1) 448 write(stdOut,*) 'Setting local fully diagonal (FD) BP dephasing model' 449 !write(stdOut,*) 'coupling=',elph%coupling 450 call set_bp_dephasing(negf, elph%coupling) 451 case(2) 452 write(stdOut,*) 'Setting local block diagonal (BD) BP dephasing model' 453 call error('NOT IMPLEMENTED! INTERRUPTED!') 454 case(3) 455 write(stdOut,*) 'Setting overlap mask (OM) block diagonal BP dephasing model' 456 call error('NOT IMPLEMENTED! INTERRUPTED!') 457 case default 458 call error("BP model is not supported") 459 end select 460 461 end subroutine negf_init_bp 462 463 !> Initialise compressed sparse row matrices 464 subroutine negf_init_csr(iAtomStart, iNeighbor, nNeighbor, img2CentCell, orb) 465 466 !> Start of orbitals for each atom 467 integer, intent(in) :: iAtomStart(:) 468 469 !> neighbours of each atom 470 integer, intent(in) :: iNeighbor(0:,:) 471 472 !> number of neighbours for each atom 473 integer, intent(in) :: nNeighbor(:) 474 475 !> mapping from image atoms to central cell 476 integer, intent(in) :: img2CentCell(:) 477 478 !> atomic orbital information 479 type(TOrbitals), intent(in) :: orb 480 481 pCsrHam => csrHam 482 pCsrOver => csrOver 483 if (allocated(csrHam%nzval)) then 484 call destroy(csrHam) 485 end if 486 call init(csrHam, iAtomStart, iNeighbor, nNeighbor, img2CentCell, orb) 487 if (allocated(csrOver%nzval)) then 488 call destroy(csrOver) 489 end if 490 call init(csrOver, csrHam) 491 492 end subroutine negf_init_csr 493 494 495 !> Destroy (module stored!) CSR matrices 496 subroutine negf_destroy() 497 498 call destruct(csrHam) 499 call destruct(csrOver) 500 call destroy_negf(negf) 501 502 write(stdOut, *) 503 write(stdOut, *) 'Release NEGF memory:' 504 !if (tIoProc) then 505 ! call writePeakInfo(6) 506 ! call writeMemInfo(6) 507 !end if 508 call writePeakInfo(stdOut) 509 call writeMemInfo(stdOut) 510 511 end subroutine negf_destroy 512 513 514 !> Initialise the structures for the libNEGF library 515 subroutine negf_init_str(denseDescr, transpar, greendens, iNeigh, nNeigh, img2CentCell) 516 517 !> Dense matrix information 518 Type(TDenseDescr), intent(in) :: denseDescr 519 520 !> transport calculation parameters 521 Type(TTranspar), intent(in) :: transpar 522 523 !> Green's function calculational parameters 524 Type(TNEGFGreenDensInfo) :: greendens 525 526 !> number of neighbours for each atom 527 Integer, intent(in) :: nNeigh(:) 528 529 !> mapping from image atoms to central cell 530 Integer, intent(in) :: img2CentCell(:) 531 532 !> neighbours of each atom 533 Integer, intent(in) :: iNeigh(0:,:) 534 535 integer, allocatable :: PL_end(:), cont_end(:), surf_end(:), cblk(:), ind(:) 536 integer, allocatable :: atomst(:), plcont(:) 537 integer, allocatable :: minv(:,:) 538 Integer :: natoms, ncont, nbl, iatc1, iatc2, iatm2 539 integer :: i, m, i1, j1, info 540 integer, allocatable :: inRegion(:) 541 542 iatm2 = transpar%idxdevice(2) 543 ncont = transpar%ncont 544 nbl = 0 545 546 if (transpar%defined) then 547 nbl = transpar%nPLs 548 else if (greendens%defined) then 549 nbl = greendens%nPLs 550 endif 551 552 if (nbl.eq.0) then 553 call error('Internal ERROR: nbl = 0 ?!') 554 end if 555 556 natoms = size(denseDescr%iatomstart) - 1 557 558 call check_pls(transpar, greendens, natoms, iNeigh, nNeigh, img2CentCell, info) 559 560 allocate(PL_end(nbl)) 561 allocate(atomst(nbl+1)) 562 allocate(plcont(nbl)) 563 allocate(cblk(ncont)) 564 allocate(cont_end(ncont)) 565 allocate(surf_end(ncont)) 566 allocate(ind(natoms+1)) 567 allocate(minv(nbl,ncont)) 568 569 ind(:) = DenseDescr%iatomstart(:) - 1 570 minv = 0 571 cblk = 0 572 573 do i = 1, ncont 574 cont_end(i) = ind(transpar%contacts(i)%idxrange(2)+1) 575 surf_end(i) = ind(transpar%contacts(i)%idxrange(1)) 576 enddo 577 578 if (transpar%defined) then 579 do i = 1, nbl-1 580 PL_end(i) = ind(transpar%PL(i+1)) 581 enddo 582 atomst(1:nbl) = transpar%PL(1:nbl) 583 PL_end(nbl) = ind(transpar%idxdevice(2)+1) 584 atomst(nbl+1) = iatm2 + 1 585 else if (greendens%defined) then 586 do i = 1, nbl-1 587 PL_end(i) = ind(greendens%PL(i+1)) 588 enddo 589 atomst(1:nbl) = greendens%PL(1:nbl) 590 PL_end(nbl) = ind(natoms+1) 591 atomst(nbl+1) = natoms + 1 592 endif 593 594 if (transpar%defined .and. ncont.gt.0) then 595 596 if(.not.transpar%tNoGeometry) then 597 598 ! For each PL finds the min atom index among the atoms in each contact 599 ! At the end the array minv(iPL,iCont) can have only one value != 0 600 ! for each contact and this is the interacting PL 601 ! NOTE: the algorithm works with the asymmetric neighbor-map of dftb+ 602 ! because atoms in contacts have larger indices than in the device 603 do m = 1, transpar%nPLs 604 ! Loop over all PL atoms 605 do i = atomst(m), atomst(m+1)-1 606 607 ! Loop over all contacts 608 do j1 = 1, ncont 609 610 iatc1 = transpar%contacts(j1)%idxrange(1) 611 iatc2 = transpar%contacts(j1)%idxrange(2) 612 613 i1 = minval(img2CentCell(iNeigh(1:nNeigh(i),i)), & 614 & mask = (img2CentCell(iNeigh(1:nNeigh(i),i)).ge.iatc1 .and. & 615 & img2CentCell(iNeigh(1:nNeigh(i),i)).le.iatc2) ) 616 617 if (i1.ge.iatc1 .and. i1.le.iatc2) then 618 minv(m,j1) = j1 619 endif 620 621 end do 622 end do 623 end do 624 625 626 do j1 = 1, ncont 627 628 if (all(minv(:,j1) == 0)) then 629 write(stdOut,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 630 write(stdOut,*) 'WARNING: contact',j1,' does not interact with any PL' 631 write(stdOut,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 632 minv(1,j1) = j1 633 end if 634 635 if (count(minv(:,j1).eq.j1).gt.1) then 636 write(stdOut,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 637 write(stdOut,*) 'ERROR: contact',j1,' interacts with more than one PL' 638 write(stdOut,*) ' check structure and increase PL size ' 639 write(stdOut,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 640 call error("") 641 end if 642 643 do m = 1, transpar%nPLs 644 if (minv(m,j1).eq.j1) then 645 cblk(j1) = m 646 end if 647 end do 648 649 end do 650 651 else 652 653 cblk=transpar%cblk 654 655 end if 656 657 write(stdOut,*) ' Structure info:' 658 write(stdOut,*) ' Number of PLs:',nbl 659 write(stdOut,*) ' PLs coupled to contacts:',cblk(1:ncont) 660 write(stdOut,*) 661 662 end if 663 664 call init_structure(negf, ncont, cont_end, surf_end, nbl, PL_end, cblk) 665 666 deallocate(PL_end) 667 deallocate(plcont) 668 deallocate(atomst) 669 deallocate(cblk) 670 deallocate(cont_end) 671 deallocate(surf_end) 672 deallocate(ind) 673 deallocate(minv) 674 675 end subroutine negf_init_str 676 677 678 !> Subroutine to check the principal layer (PL) definitions 679 subroutine check_pls(transPar, greenDens, nAtoms, iNeigh, nNeigh, img2CentCell, info) 680 681 !> transport calculation parameters 682 type(TTranspar), intent(in) :: transPar 683 684 !> Green's function calculational parameters 685 type(TNEGFGreenDensInfo) :: greenDens 686 687 !> Number of atoms in the system 688 integer, intent(in) :: nAtoms 689 690 !> neighbours of each atom 691 integer, intent(in) :: iNeigh(0:,:) 692 693 !> number of neighbours for each atom 694 integer, intent(in) :: nNeigh(:) 695 696 !> mapping from image atoms to central cell 697 integer, intent(in) :: img2CentCell(:) 698 699 !> error message 700 integer, intent(out) :: info 701 702 integer :: nbl, iatm1, iatm2, iats, iate 703 integer :: mm, nn, ii, kk 704 integer, allocatable :: atomst(:) 705 706 ! The contacts have been already checked 707 ! Here checks the PL definition and contact/device interactions 708 709 iatm1 = transpar%idxdevice(1) 710 iatm2 = transpar%idxdevice(2) 711 712 nbl = 0 713 714 if (transpar%defined) then 715 nbl = transpar%nPLs 716 else if (greendens%defined) then 717 nbl = greendens%nPLs 718 endif 719 720 if (nbl.eq.0) then 721 call error('Internal ERROR: nbl = 0 ?!') 722 end if 723 724 allocate(atomst(nbl+1)) 725 726 if (transpar%defined) then 727 atomst(1:nbl) = transpar%PL(1:nbl) 728 atomst(nbl+1) = iatm2 + 1 729 else if (greendens%defined) then 730 atomst(1:nbl) = greendens%PL(1:nbl) 731 atomst(nbl+1) = natoms + 1 732 endif 733 734 info = 0 735 do mm = 1, nbl-1 736 do nn = mm+1, nbl 737 iats = atomst(nn) 738 iate = atomst(nn+1)-1 739 do ii = atomst(mm), atomst(mm+1)-1 740 kk = maxval( img2CentCell(iNeigh(1:nNeigh(ii),ii)), & 741 mask = (img2CentCell(iNeigh(1:nNeigh(ii),ii)).ge.iats .and. & 742 img2CentCell(iNeigh(1:nNeigh(ii),ii)).le.iate) ) 743 end do 744 if (nn .gt. mm+1 .and. kk .ge. iats .and. kk .le. iate) then 745 write(stdOut,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 746 write(stdOut,*) 'WARNING: PL ',mm,' interacts with PL',nn 747 write(stdOut,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 748 info = mm 749 end if 750 end do 751 end do 752 753 deallocate(atomst) 754 755 end subroutine check_pls 756 757 758 !> Interface subroutine to call Density matrix computation 759 subroutine negf_density(miter,spin,nkpoint,HH,SS,mu,DensMat,EnMat) 760 761 !> SCC step (used in SGF) 762 integer, intent (in) :: miter 763 764 !> spin component (SGF) 765 integer, intent (in) :: spin 766 767 !> nk point (used in SGF) 768 integer, intent (in) :: nkpoint 769 770 !> Hamiltonian 771 type(z_CSR), pointer, intent(in) :: HH 772 773 !> Overlap 774 type(z_CSR), pointer, intent(in) :: SS 775 776 !> chemical potential 777 real(dp), intent(in) :: mu(:) 778 779 !> Density matrix 780 type(z_CSR), pointer, intent(in), optional :: DensMat 781 782 !> Energy weighted DM 783 type(z_CSR), pointer, intent(in), optional :: EnMat 784 785 type(lnParams) :: params 786 integer :: nn 787 788 call get_params(negf, params) 789 790 params%kpoint = nkpoint 791 params%spin = spin 792 params%DorE='N' 793 nn=size(mu,1) 794 params%mu(1:nn) = mu(1:nn) 795 796 if(present(DensMat)) then 797 params%DorE = 'D' 798 call set_params(negf,params) 799 call pass_DM(negf,rho=DensMat) 800 endif 801 if(present(EnMat)) then 802 params%DorE = 'E' 803 call set_params(negf,params) 804 call pass_DM(negf,rhoE=EnMat) 805 endif 806 if (present(DensMat).and.present(EnMat)) then 807 params%DorE = 'B' 808 call set_params(negf,params) 809 call error('UNSUPPORTED CASE in negf_density') 810 endif 811 812 if (params%DorE.eq.'N') then 813 return 814 end if 815 816 call pass_HS(negf,HH,SS) 817 818 call compute_density_dft(negf) 819 820 call destroy_matrices(negf) 821 822 end subroutine negf_density 823 824 825 !> Interface subroutine to call ldos computation 826 subroutine negf_ldos(HH, SS, spin, kpoint, wght, ledos) 827 828 !> hamiltonian in CSR format 829 type(z_CSR), pointer, intent(in) :: HH 830 831 !> overlap in CSR format 832 type(z_CSR), pointer, intent(in) :: SS 833 834 !> spin index 835 integer, intent(in) :: spin 836 837 !> kpoint index 838 integer, intent(in) :: kpoint 839 840 !> kpoint weight 841 real(dp), intent(in) :: wght 842 843 !> local DOS 844 real(dp), dimension(:,:), pointer :: ledos 845 846 type(lnParams) :: params 847 848 call get_params(negf, params) 849 850 params%spin = spin 851 params%kpoint = kpoint 852 params%wght = wght 853 854 call pass_HS(negf,HH,SS) 855 856 call compute_ldos(negf) 857 858 call destroy_matrices(negf) 859 860 call associate_ldos(pNegf, ledos) 861 862 end subroutine negf_ldos 863 864 865 !> Debug routine to dump H and S as a file in Matlab format 866 subroutine negf_dumpHS(HH,SS) 867 868 !> hamiltonian in CSR format 869 type(z_CSR), intent(in) :: HH 870 871 !> Overlap in CSR format 872 type(z_CSR), intent(in) :: SS 873 874 integer :: fdUnit 875 876 write(stdOut, *) 'Dumping H and S in files...' 877 878 open(newunit=fdUnit, file='HH.dat') 879 write(fdUnit, *) '% Size =',HH%nrow, HH%ncol 880 write(fdUnit, *) '% Nonzeros =',HH%nnz 881 write(fdUnit, *) '% ' 882 write(fdUnit, *) 'zzz = [' 883 call printcsr(fdUnit, HH) 884 write(fdUnit, *) ']' 885 close(fdUnit) 886 887 open(newunit=fdUnit, file='SS.dat') 888 write(fdUnit, *) '% Size =',SS%nrow, SS%ncol 889 write(fdUnit, *) '% Nonzeros =',SS%nnz 890 write(fdUnit, *) '% ' 891 write(fdUnit, *) 'zzz = [' 892 call printcsr(fdUnit, SS) 893 write(fdUnit, *) ']' 894 close(fdUnit) 895 896 end subroutine negf_dumpHS 897 898 899 !> Routines to setup orthogonalised H and S have been moved here 900 subroutine prepare_HS(H_dev,S_dev,HH,SS) 901 902 !> hamiltonian in dense format 903 real(dp), dimension(:,:) :: H_dev 904 905 !> overlap in dense format 906 real(dp), dimension(:,:) :: S_dev 907 908 !> hamiltonian in CSR format 909 type(z_CSR), intent(inout) :: HH 910 911 !> overlap in CSR format 912 type(z_CSR), intent(inout) :: SS 913 914 if (negf%tOrthonormal) then 915 write(stdOut, "(' Lowdin orthogonalization for the whole system ')") 916 call Orthogonalization(H_dev, S_dev) 917 end if 918 919 if (negf%tOrthonormalDevice) then 920 write(stdOut, "(' Lowdin orthogonalization for device-only')") 921 call Orthogonalization_dev(H_dev, S_dev) 922 end if 923 924 if (negf%tOrthonormal.or.negf%tOrthonormalDevice) then 925 call MakeHHSS(H_dev,S_dev,HH,SS) 926 end if 927 928 !if(negf%tManyBody) call MakeHS_dev 929 930 end subroutine 931 932 933 !> Interface subroutine to call calculation of currents 934 subroutine negf_current(HH, SS, spin, kpoint, wght, tunn, curr, ledos, currents) 935 936 !> hamiltonian 937 type(z_CSR), pointer, intent(in) :: HH 938 939 !> overlap 940 type(z_CSR), pointer, intent(in) :: SS 941 942 !> spin index 943 integer, intent(in) :: spin 944 945 !> kpoint index 946 integer, intent(in) :: kpoint 947 948 !> kpoint weight 949 real(dp), intent(in) :: wght 950 951 !> Tunneling amplitudes 952 real(dp), dimension(:,:), pointer :: tunn 953 954 !> current magnitudes 955 real(dp), dimension(:,:), pointer :: curr 956 957 !> local density of states 958 real(dp), dimension(:,:), pointer :: ledos 959 960 !> current directions 961 real(dp), dimension(:), pointer :: currents 962 963 type(lnParams) :: params 964 965 call get_params(negf, params) 966 967 params%spin = spin 968 params%kpoint = kpoint 969 params%wght = wght 970 971 call set_params(negf, params) 972 973 call pass_HS(negf,HH,SS) 974 975 call compute_current(negf) 976 977 ! Associate internal negf arrays to local pointers 978 call associate_ldos(pNegf, ledos) 979 call associate_transmission(pNegf, tunn) 980 call associate_current(pNegf, curr) 981 982 call associate_lead_currents(pNegf, currents) 983 if (.not.associated(currents)) then 984 call error('Internal error: currVec not associated') 985 end if 986 987 end subroutine negf_current 988 989 990 991 !> Calculates density matrix with Green's functions 992#:if WITH_MPI 993 subroutine calcdensity_green(iSCCIter, mpicomm, groupKS, ham, over, iNeighbor, nNeighbor,& 994 & iAtomStart, iPair, img2CentCell, iCellVec, cellVec, orb, kPoints, kWeights, mu, rho, Eband,& 995 & Ef, E0, TS) 996 997 !> MPI communicator 998 type(mpifx_comm), intent(in) :: mpicomm 999#:else 1000 subroutine calcdensity_green(iSCCIter, groupKS, ham, over, iNeighbor, nNeighbor,& 1001 & iAtomStart, iPair, img2CentCell, iCellVec, cellVec, orb, kPoints, kWeights, mu, rho, Eband,& 1002 & Ef, E0, TS) 1003#:endif 1004 1005 !> SCC iteration 1006 integer, intent(in) :: iSCCIter 1007 1008 !> kpoint and spin descriptor 1009 integer, intent(in) :: groupKS(:,:) 1010 1011 !> hamiltonian matrix 1012 real(dp), intent(in) :: ham(:,:) 1013 1014 !> overlap matrix 1015 real(dp), intent(in) :: over(:) 1016 1017 !> neighbours of atoms 1018 integer, intent(in) :: iNeighbor(0:,:) 1019 1020 !> number of neighbours 1021 integer, intent(in) :: nNeighbor(:) 1022 1023 !> dense indexing for orbitals 1024 integer, intent(in) :: iAtomStart(:) 1025 1026 !> sparse indexing for orbitals 1027 integer, intent(in) :: iPair(0:,:) 1028 1029 !> map from image to central cell atoms 1030 integer, intent(in) :: img2CentCell(:) 1031 1032 !> index for unit cell an atom is associated with 1033 integer, intent(in) :: iCellVec(:) 1034 1035 !> Vectors to unit cells 1036 real(dp), intent(in) :: cellVec(:,:) 1037 1038 !> atomic orbital information 1039 type(TOrbitals), intent(in) :: orb 1040 1041 !> k-points 1042 real(dp), intent(in) :: kPoints(:,:) 1043 1044 !> k-point weights 1045 real(dp), intent(in) :: kWeights(:) 1046 1047 !> chemical potentials of reservoirs 1048 real(dp), intent(in) :: mu(:,:) 1049 1050 !> density matrix 1051 real(dp), intent(out) :: rho(:,:) 1052 1053 !> band energy 1054 real(dp), intent(out) :: Eband(:) 1055 1056 !> Fermi energy 1057 real(dp), intent(out) :: Ef(:) 1058 1059 !> zero temperature (extrapolated) electronic energy 1060 real(dp), intent(out) :: E0(:) 1061 1062 !> Electron entropy 1063 real(dp), intent(out) :: TS(:) 1064 1065 integer :: nSpin, nKS, iK, iS, iKS 1066 type(z_CSR), target :: csrDens 1067 type(z_CSR), pointer :: pCsrDens 1068 1069 pCsrDens => csrDens 1070 1071#:if WITH_MPI 1072 call negf_mpi_init(mpicomm) 1073#:endif 1074 !Decide what to do with surface GFs. 1075 !sets readOldSGF: if it is 0 or 1 it is left so 1076 if (negf%readOldDM_SGFs.eq.COMPSAVE_SGF) then 1077 if(iSCCIter.eq.1) then 1078 call set_readOldDMsgf(negf, COMPSAVE_SGF) ! compute and save SGF on files 1079 else 1080 call set_readOldDMsgf(negf, READ_SGF) ! read from files 1081 endif 1082 endif 1083 ! We need this now for different fermi levels in colinear spin 1084 ! Note: the spin polirized does not work with 1085 ! built-int potentials (the unpolarized does) in the poisson 1086 nKS = size(groupKS, dim=2) 1087 nSpin = size(ham, dim=2) 1088 rho = 0.0_dp 1089 1090 write(stdOut, *) 1091 write(stdOut, '(80("="))') 1092 write(stdOut, *) ' COMPUTING DENSITY MATRIX ' 1093 write(stdOut, '(80("="))') 1094 1095 1096 do iKS = 1, nKS 1097 iK = groupKS(1, iKS) 1098 iS = groupKS(2, iKS) 1099 1100 write(stdOut,*) 'k-point',iK,'Spin',iS 1101 1102 call foldToCSR(csrHam, ham(:,iS), kPoints(:,iK), iAtomStart, iPair, iNeighbor, nNeighbor,& 1103 & img2CentCell, iCellVec, cellVec, orb) 1104 call foldToCSR(csrOver, over, kPoints(:,ik), iAtomStart, iPair, iNeighbor, nNeighbor,& 1105 & img2CentCell, iCellVec, cellVec, orb) 1106 1107 call negf_density(iSCCIter, iS, iKS, pCsrHam, pCsrOver, mu(:,iS), DensMat=pCsrDens) 1108 1109 ! NOTE: 1110 ! unfold adds up to rho the csrDens(k) contribution 1111 call unfoldFromCSR(rho(:,iS), csrDens, kPoints(:,iK), kWeights(iK), iAtomStart, iPair,& 1112 & iNeighbor, nNeighbor, img2CentCell, iCellVec, cellVec, orb) 1113 1114 call destruct(csrDens) 1115 1116 ! Set some fake energies: 1117 Eband(iS) = 0.0_dp 1118 Ef(iS) = mu(1,iS) 1119 TS(iS) = 0.0_dp 1120 E0(iS) = 0.0_dp 1121 1122 end do 1123 1124#:if WITH_MPI 1125 do iS = 1, nSpin 1126 ! In place all-reduce of the density matrix 1127 call mpifx_allreduceip(mpicomm, rho(:,iS), MPI_SUM) 1128 end do 1129#:endif 1130 1131 ! Now SGFs can be read unless not stored 1132 if (negf%readOldDM_SGFs.ne.COMP_SGF) then 1133 call set_readOldDMsgf(negf, READ_SGF) ! read from files 1134 end if 1135 1136 write(stdOut,'(80("="))') 1137 write(stdOut,*) 1138 1139 end subroutine calcdensity_green 1140 1141 !> Calculates energy-weighted density matrix with Green's functions 1142#:if WITH_MPI 1143 subroutine calcEdensity_green(iSCCIter, mpicomm, groupKS, ham, over, iNeighbor, nNeighbor,& 1144 & iAtomStart, iPair, img2CentCell, iCellVec, cellVec, orb, kPoints, kWeights, mu, rhoE) 1145 1146 !> MPI communicator 1147 type(mpifx_comm), intent(in) :: mpicomm 1148#:else 1149 subroutine calcEdensity_green(iSCCIter, groupKS, ham, over, iNeighbor, nNeighbor,& 1150 & iAtomStart, iPair, img2CentCell, iCellVec, cellVec, orb, kPoints, kWeights, mu, rhoE) 1151#:endif 1152 1153 !> SCC iteration 1154 integer, intent(in) :: iSCCIter 1155 1156 !> kpoint and spin descriptor 1157 integer, intent(in) :: groupKS(:,:) 1158 1159 !> hamiltonian matrix 1160 real(dp), intent(in) :: ham(:,:) 1161 1162 !> overlap matrix 1163 real(dp), intent(in) :: over(:) 1164 1165 !> neighbours of atoms 1166 integer, intent(in) :: iNeighbor(0:,:) 1167 1168 !> number of neighbours 1169 integer, intent(in) :: nNeighbor(:) 1170 1171 !> dense indexing for orbitals 1172 integer, intent(in) :: iAtomStart(:) 1173 1174 !> sparse indexing for orbitals 1175 integer, intent(in) :: iPair(0:,:) 1176 1177 !> map from image to central cell atoms 1178 integer, intent(in) :: img2CentCell(:) 1179 1180 !> index for unit cell an atom is associated with 1181 integer, intent(in) :: iCellVec(:) 1182 1183 !> Vectors to unit cells 1184 real(dp), intent(in) :: cellVec(:,:) 1185 1186 !> atomic orbital information Needs only orb%nOrbAtom, orb%mOrb 1187 type(TOrbitals), intent(in) :: orb 1188 1189 !> k-points 1190 real(dp), intent(in) :: kPoints(:,:) 1191 1192 !> k-point weights 1193 real(dp), intent(in) :: kWeights(:) 1194 1195 !> chemical potentials 1196 real(dp), intent(in) :: mu(:,:) 1197 1198 !> Energy weighted density matrix 1199 real(dp), intent(out) :: rhoE(:) 1200 1201 integer :: nSpin, nKS, iK, iS, iKS 1202 type(z_CSR), target :: csrEDens 1203 type(z_CSR), pointer :: pCsrEDens 1204 1205 pCsrEDens => csrEDens 1206 1207#:if WITH_MPI 1208 call negf_mpi_init(mpicomm) 1209#:endif 1210 !Decide what to do with surface GFs. 1211 !sets readOldSGF: if it is 0 or 1 it is left so 1212 if (negf%readOldDM_SGFs.eq.COMPSAVE_SGF) then 1213 if(iSCCIter.eq.1) then 1214 call set_readOldDMsgf(negf, COMPSAVE_SGF) ! compute and save SGF on files 1215 else 1216 call set_readOldDMsgf(negf, READ_SGF) ! read from files 1217 endif 1218 endif 1219 1220 ! We need this now for different fermi levels in colinear spin 1221 ! Note: the spin polirized does not work with 1222 ! built-int potentials (the unpolarized does) in the poisson 1223 ! I do not set the fermi because it seems that in libnegf it is 1224 ! not really needed 1225 1226 nKS = size(groupKS, dim=2) 1227 nSpin = size(ham, dim=2) 1228 rhoE = 0.0_dp 1229 1230 write(stdOut, *) 1231 write(stdOut, '(80("="))') 1232 write(stdOut, *) ' COMPUTING E-WEIGHTED DENSITY MATRIX ' 1233 write(stdOut, '(80("="))') 1234 1235 do iKS = 1, nKS 1236 iK = groupKS(1, iKS) 1237 iS = groupKS(2, iKS) 1238 1239 write(stdOut,*) 'k-point',iK,'Spin',iS 1240 1241 call foldToCSR(csrHam, ham(:,iS), kPoints(:,iK), iAtomStart, iPair, iNeighbor, nNeighbor,& 1242 & img2CentCell, iCellVec, cellVec, orb) 1243 call foldToCSR(csrOver, over, kPoints(:,ik), iAtomStart, iPair, iNeighbor, nNeighbor,& 1244 & img2CentCell, iCellVec, cellVec, orb) 1245 1246 call negf_density(iSCCIter, iS, iKS, pCsrHam, pCsrOver, mu(:,iS), EnMat=pCsrEDens) 1247 1248 ! NOTE: 1249 ! unfold adds up to rhoEPrim the csrEDens(k) contribution 1250 ! 1251 call unfoldFromCSR(rhoE, csrEDens, kPoints(:,iK), kWeights(iK), iAtomStart, iPair, iNeighbor,& 1252 & nNeighbor, img2CentCell, iCellVec, cellVec, orb) 1253 1254 call destruct(csrEDens) 1255 1256 end do 1257 1258 ! In place all-reduce of the density matrix 1259#:if WITH_MPI 1260 call mpifx_allreduceip(mpicomm, rhoE, MPI_SUM) 1261#:endif 1262 1263 ! Now SGFs can be read unless not stored 1264 if (negf%readOldDM_SGFs.ne.COMP_SGF) then 1265 call set_readOldDMsgf(negf, READ_SGF) ! read from files 1266 end if 1267 1268 write(stdOut,'(80("="))') 1269 write(stdOut,*) 1270 1271 end subroutine calcEdensity_green 1272 1273 1274 !> Calculate the current and optionally density of states 1275#:if WITH_MPI 1276 subroutine calc_current(mpicomm, groupKS, ham, over, iNeighbor, nNeighbor, iAtomStart, iPair,& 1277 & img2CentCell, iCellVec, cellVec, orb, kPoints, kWeights, tunnMat, currMat, ldosMat,& 1278 & currLead, writeTunn, tWriteLDOS, regionLabelLDOS, mu) 1279 1280 !> MPI communicator 1281 type(mpifx_comm), intent(in) :: mpicomm 1282#:else 1283 subroutine calc_current(groupKS, ham, over, iNeighbor, nNeighbor, iAtomStart, iPair,& 1284 & img2CentCell, iCellVec, cellVec, orb, kPoints, kWeights, tunnMat, currMat, ldosMat,& 1285 & currLead, writeTunn, tWriteLDOS, regionLabelLDOS, mu) 1286#:endif 1287 1288 !> kpoint and spin descriptor 1289 integer, intent(in) :: groupKS(:,:) 1290 1291 !> hamiltonian matrix 1292 real(dp), intent(in) :: ham(:,:) 1293 1294 !> overlap matrix 1295 real(dp), intent(in) :: over(:) 1296 1297 !> neighbours of atoms 1298 integer, intent(in) :: iNeighbor(0:,:) 1299 1300 !> number of neighbours 1301 integer, intent(in) :: nNeighbor(:) 1302 1303 !> dense indexing for orbitals 1304 integer, intent(in) :: iAtomStart(:) 1305 1306 !> sparse indexing for orbitals 1307 integer, intent(in) :: iPair(0:,:) 1308 1309 !> map from image to central cell atoms 1310 integer, intent(in) :: img2CentCell(:) 1311 1312 !> index for unit cell an atom is associated with 1313 integer, intent(in) :: iCellVec(:) 1314 1315 !> Vectors to unit cells 1316 real(dp), intent(in) :: cellVec(:,:) 1317 1318 !> atomic orbital information Needs only orb%nOrbAtom, orb%mOrb 1319 type(TOrbitals), intent(in) :: orb 1320 1321 !> k-points 1322 real(dp), intent(in) :: kPoints(:,:) 1323 1324 !> k-point weights 1325 real(dp), intent(in) :: kWeights(:) 1326 1327 !> matrix of tunnelling amplitudes at each energy from contacts 1328 real(dp), allocatable, intent(inout) :: tunnMat(:,:) 1329 1330 !> matrix of current at each energy from currents 1331 real(dp), allocatable, intent(inout) :: currMat(:,:) 1332 1333 !> density of states for each energy and region of projection 1334 real(dp), allocatable, intent(inout) :: ldosMat(:,:) 1335 1336 !> current into/out of contacts 1337 real(dp), allocatable, intent(inout) :: currLead(:) 1338 1339 !> should tunneling data be written 1340 logical, intent(in) :: writeTunn 1341 1342 !> should DOS data be written 1343 logical, intent(in) :: tWriteLDOS 1344 1345 !> labels for DOS projected regions 1346 character(lc), allocatable, intent(in) :: regionLabelLDOS(:) 1347 1348 !> We need this now for different fermi levels in colinear spin 1349 !> Note: the spin polarised does not work with 1350 !> built-int potentials (the unpolarized does) in the poisson 1351 !> I do not set the fermi because it seems that in libnegf it is 1352 !> not really needed 1353 real(dp), intent(in) :: mu(:,:) 1354 1355 real(dp), allocatable :: tunnSKRes(:,:,:), currSKRes(:,:,:), ldosSKRes(:,:,:) 1356 real(dp), pointer :: tunnPMat(:,:)=>null() 1357 real(dp), pointer :: currPMat(:,:)=>null() 1358 real(dp), pointer :: ldosPMat(:,:)=>null() 1359 real(dp), pointer :: currPVec(:)=>null() 1360 integer :: iKS, iK, iS, nKS, ii, err, ncont, readSGFbkup 1361 type(unit) :: unitOfEnergy ! Set the units of H 1362 type(unit) :: unitOfCurrent ! Set desired units for Jel 1363 type(lnParams) :: params 1364 1365 integer :: i, j, k, NumStates, icont 1366 real(dp), dimension(:,:), allocatable :: H_all, S_all 1367 character(:), allocatable :: filename 1368 1369#:if WITH_MPI 1370 call negf_mpi_init(mpicomm) 1371#:endif 1372 call get_params(negf, params) 1373 1374 unitOfEnergy%name = "H" 1375 unitOfCurrent%name = "A" 1376 1377 nKS = size(groupKS, dim=2) 1378 ncont = size(mu,1) 1379 1380 if (params%verbose.gt.30) then 1381 write(stdOut, *) 1382 write(stdOut, '(80("="))') 1383 write(stdOut, *) ' COMPUTATION OF CURRENT ' 1384 write(stdOut, '(80("="))') 1385 write(stdOut, *) 1386 end if 1387 1388 do iKS = 1, nKS 1389 iK = groupKS(1, iKS) 1390 iS = groupKS(2, iKS) 1391 1392 write(stdOut,*) 'Spin',iS,'k-point',iK,'k-weight',kWeights(iK) 1393 1394 params%mu(1:ncont) = mu(1:ncont,iS) 1395 1396 call set_params(negf, params) 1397 1398 if (negf%NumStates.eq.0) then 1399 negf%NumStates=csrHam%ncol 1400 end if 1401 1402 !*** ORTHOGONALIZATIONS *** 1403 ! THIS MAKES SENSE ONLY FOR A REAL MATRICES, i.e. k==0 && collinear spin 1404 if (all(kPoints(:,iK) .eq. 0.0_dp) .and. & 1405 (negf%tOrthonormal .or. negf%tOrthonormalDevice)) then 1406 1407 NumStates = negf%NumStates 1408 1409 if (.not.allocated(H_all)) then 1410 allocate(H_all(NumStates,NumStates)) 1411 end if 1412 if (.not.allocated(S_all)) then 1413 allocate(S_all(NumStates,NumStates)) 1414 end if 1415 1416 call unpackHS(H_all, ham(:,iS), iNeighbor, nNeighbor, iAtomStart, iPair, img2CentCell) 1417 call blockSymmetrizeHS(H_all, iAtomStart) 1418 1419 call unpackHS(S_all, over, iNeighbor, nNeighbor, iAtomStart, iPair, img2CentCell) 1420 call blockSymmetrizeHS(S_all, iAtomStart) 1421 1422 call prepare_HS(H_all,S_all,csrHam,csrOver) 1423 1424 else 1425 1426 call foldToCSR(csrHam, ham(:,iS), kPoints(:,iK), iAtomStart, iPair, iNeighbor, nNeighbor,& 1427 & img2CentCell, iCellVec, cellVec, orb) 1428 1429 call foldToCSR(csrOver, over, kPoints(:,ik), iAtomStart, iPair, iNeighbor, nNeighbor,& 1430 & img2CentCell, iCellVec, cellVec, orb) 1431 1432 end if 1433 1434 call negf_current(pCsrHam, pCsrOver, iS, iK, kWeights(iK), tunnPMat, currPMat, ldosPMat,& 1435 & currPVec) 1436 1437 if(.not.allocated(currLead)) then 1438 allocate(currLead(size(currPVec)), stat=err) 1439 if (err /= 0) then 1440 call error('Allocation error (currTot)') 1441 end if 1442 currLead = 0.0_dp 1443 endif 1444 currLead = currLead + currPVec 1445 1446 !GUIDE: tunnPMat libNEGF output stores Transmission, T(iE, i->j) 1447 ! tunnPMat is MPI distributed on energy points (0.0 on other nodes) 1448 ! tunnMat MPI gather partial results and accumulate k-summation 1449 ! currPMat stores contact current I_i(iE) 1450 ! tunnSKRes stores tunneling for all k-points and spin: T(iE, i->j, iSK) 1451 1452#:if WITH_MPI 1453 call add_partial_results(mpicomm, tunnPMat, tunnMat, tunnSKRes, iKS, nKS) 1454 call add_partial_results(mpicomm, currPMat, currMat, currSKRes, iKS, nKS) 1455 call add_partial_results(mpicomm, ldosPMat, ldosMat, ldosSKRes, iKS, nKS) 1456#:else 1457 call add_partial_results(tunnPMat, tunnMat, tunnSKRes, iKS, nKS) 1458 call add_partial_results(currPMat, currMat, currSKRes, iKS, nKS) 1459 call add_partial_results(ldosPMat, ldosMat, ldosSKRes, iKS, nKS) 1460#:endif 1461 1462 end do 1463 1464#:if WITH_MPI 1465 call mpifx_allreduceip(mpicomm, currLead, MPI_SUM) 1466 call mpifx_barrier(mpicomm) 1467#:endif 1468 1469 ! converts from internal atomic units into amperes 1470 currLead = currLead * convertCurrent(unitOfEnergy, unitOfCurrent) 1471 1472 do ii = 1, size(currLead) 1473 write(stdOut, *) 1474 write(stdOut, '(1x,a,i3,i3,a,ES14.5,a,a)') ' contacts: ',params%ni(ii),params%nf(ii),& 1475 & ' current: ', currLead(ii),' ',unitOfCurrent%name 1476 enddo 1477 1478 ! Write Total transmission, T(E), on a separate file (optional) 1479 if (allocated(tunnMat)) then 1480 filename = 'transmission' 1481 if (tIOProc .and. writeTunn) then 1482 call write_file(negf, tunnMat, tunnSKRes, filename, groupKS, kpoints, kWeights) 1483 end if 1484 if (allocated(tunnSKRes)) then 1485 deallocate(tunnSKRes) 1486 end if 1487 else 1488 ! needed to avoid some segfault 1489 if (allocated(tunnMat)) then 1490 deallocate(tunnMat) 1491 end if 1492 allocate(tunnMat(0,0)) 1493 end if 1494 1495 ! Write Total lead current, I_i(E), on a separate file (optional) 1496 if (allocated(currMat)) then 1497 filename = 'current' 1498 if (tIOProc .and. writeTunn) then 1499 call write_file(negf, currMat, currSKRes, filename, groupKS, kpoints, kWeights) 1500 end if 1501 if (allocated(currSKRes)) then 1502 deallocate(currSKRes) 1503 end if 1504 else 1505 ! needed to avoid some segfault 1506 if (allocated(currMat)) then 1507 deallocate(currMat) 1508 end if 1509 allocate(currMat(0,0)) 1510 endif 1511 1512 if (allocated(ldosMat)) then 1513 ! Write Total localDOS on a separate file (optional) 1514 if (tIoProc .and. tWriteLDOS) then 1515 @:ASSERT(allocated(regionLabelLDOS)) 1516 call write_files(negf, ldosMat, ldosSKRes, groupKS, kpoints, kWeights, regionLabelLDOS) 1517 if (allocated(ldosSKRes)) then 1518 deallocate(ldosSKRes) 1519 end if 1520 else 1521 ! needed to avoid some segfault 1522 if (allocated(ldosMat)) then 1523 deallocate(ldosMat) 1524 end if 1525 allocate(ldosMat(0,0)) 1526 end if 1527 end if 1528 1529 end subroutine calc_current 1530 1531 !> utility to allocate and sum partial results from different channels 1532#:if WITH_MPI 1533 subroutine add_partial_results(mpicomm, pMat, matTot, matSKRes, iK, nK) 1534 1535 !> MPI communicator 1536 type(mpifx_comm), intent(in) :: mpicomm 1537#:else 1538 subroutine add_partial_results(pMat, matTot, matSKRes, iK, nK) 1539#:endif 1540 1541 !> pointer to matrix of data 1542 real(dp), intent(in), pointer :: pMat(:,:) 1543 1544 !> sum total 1545 real(dp), allocatable, intent(inout) :: matTot(:,:) 1546 1547 !> k-resolved sum 1548 real(dp), allocatable, intent(inout) :: matSKRes(:,:,:) 1549 1550 !> particular k-point 1551 integer, intent(in) :: iK 1552 1553 !> number of k-points 1554 integer, intent(in) :: nK 1555 1556 real(dp), allocatable :: tmpMat(:,:) 1557 integer :: err 1558 1559 if (associated(pMat)) then 1560#:if WITH_MPI 1561 allocate(tmpMat(size(pMat,1), size(pMat,2)), stat=err) 1562 1563 if (err /= 0) then 1564 call error('Allocation error (tmpMat)') 1565 end if 1566 1567 tmpMat = pMat 1568 call mpifx_allreduceip(mpicomm, tmpMat, MPI_SUM) 1569#:endif 1570 if(.not.allocated(matTot)) then 1571 allocate(matTot(size(pMat,1), size(pMat,2)), stat=err) 1572 1573 if (err /= 0) then 1574 call error('Allocation error (tunnTot)') 1575 end if 1576 1577 matTot = 0.0_dp 1578 end if 1579#:if WITH_MPI 1580 matTot = matTot + tmpMat 1581#:else 1582 matTot = matTot + pMat 1583#:endif 1584 1585 if (nK > 1) then 1586 if (.not.allocated(matSKRes)) then 1587 allocate(matSKRes(size(pMat,1), size(pMat,2), nK), stat=err) 1588 1589 if (err/=0) then 1590 call error('Allocation error (tunnSKRes)') 1591 end if 1592 1593 matSKRes = 0.0_dp 1594 endif 1595#:if WITH_MPI 1596 matSKRes(:,:,iK) = tmpMat(:,:) 1597#:else 1598 matSKRes(:,:,iK) = pMat(:,:) 1599#:endif 1600 end if 1601 1602#:if WITH_MPI 1603 deallocate(tmpMat) 1604#:endif 1605 1606 end if 1607 1608 end subroutine add_partial_results 1609 1610 !> utility to write tunneling files 1611 subroutine write_file(negf, matTot, matSKRes, filename, groupKS, kpoints, kWeights) 1612 1613 !> Contains input data, runtime quantities and output data 1614 type(TNegf) :: negf 1615 1616 !> results to print if allocated 1617 real(dp), intent(in), allocatable :: matTot(:,:) 1618 1619 !> k- and spin-resolved quantities, if allocated 1620 real(dp), intent(in), allocatable :: matSKRes(:,:,:) 1621 1622 !> file to print out to 1623 character(*), intent(in) :: filename 1624 1625 !> local k-points and spins on this processor 1626 integer, intent(in) :: groupKS(:,:) 1627 1628 !> k-points 1629 real(dp), intent(in) :: kPoints(:,:) 1630 1631 !> Weights for k-points 1632 real(dp), intent(in) :: kWeights(:) 1633 1634 integer :: ii, jj, nKS, iKS, nK, nS, iK, iS, fdUnit 1635 type(lnParams) :: params 1636 1637 call get_params(negf, params) 1638 1639 nKS = size(groupKS, dim=2) 1640 nK = size(kpoints, dim=2) 1641 nS = nKS/nK 1642 1643 open(newunit=fdUnit, file=trim(filename)//'.dat') 1644 do ii=1,size(matTot, dim=1) 1645 write(fdUnit,'(F20.6)',ADVANCE='NO') (params%Emin+(ii-1)*params%Estep) * Hartree__eV 1646 do jj=1,size(matTot, dim=2) 1647 write(fdUnit,'(ES20.8)',ADVANCE='NO') matTot(ii,jj) 1648 enddo 1649 write(fdUnit,*) 1650 enddo 1651 close(fdUnit) 1652 1653 if (nKS > 1) then 1654 1655 open(newunit=fdUnit, file=trim(filename)//'_kpoints.dat') 1656 write(fdUnit,*)'# NKpoints = ', nK 1657 write(fdUnit,*)'# NSpin = ', nS 1658 write(fdUnit,*)'# Energy [eV], <spin k1 k2 k3 weight> ' 1659 write(fdUnit,'(A1)', ADVANCE='NO') '# ' 1660 1661 do iKS = 1, nKS 1662 iK = groupKS(1,iKS) 1663 iS = groupKS(2,iKS) 1664 write(fdUnit,'(i5.2)', ADVANCE='NO') iS 1665 write(fdUnit,'(es15.5, es15.5, es15.5, es15.5)', ADVANCE='NO') kpoints(:,iK), kWeights(iK) 1666 end do 1667 write(fdUnit,*) 1668 1669 if (allocated(matSKRes)) then 1670 do ii = 1, size(matSKRes(:,:,:), dim=1) 1671 write(fdUnit,'(f20.6)',ADVANCE='NO') (params%Emin+(ii-1)*params%Estep) * Hartree__eV 1672 do jj = 1, size(matSKRes(:,:,:), dim=2) 1673 do iKS = 1,nKS 1674 write(fdUnit,'(es20.8)',ADVANCE='NO') matSKRes(ii,jj, iKS) 1675 enddo 1676 write(fdUnit,*) 1677 enddo 1678 enddo 1679 end if 1680 close(fdUnit) 1681 1682 end if 1683 1684 end subroutine write_file 1685 1686 1687 !> utility to write tunneling/ldos files with names from labels 1688 subroutine write_files(negf, matTot, matSKRes, groupKS, kpoints, kWeights, regionLabels) 1689 1690 !> Contains input data, runtime quantities and output data 1691 type(TNegf) :: negf 1692 1693 !> results to print if allocated 1694 real(dp), intent(in) :: matTot(:,:) 1695 1696 !> k- and spin-resolved quantities, if allocated 1697 real(dp), intent(in), allocatable :: matSKRes(:,:,:) 1698 1699 !> local k-points and spins on this processor 1700 integer, intent(in) :: groupKS(:,:) 1701 1702 !> k-points 1703 real(dp), intent(in) :: kPoints(:,:) 1704 1705 !> Weights for k-points 1706 real(dp), intent(in) :: kWeights(:) 1707 1708 !> Labels for the separate regions 1709 character(lc), intent(in) :: regionLabels(:) 1710 1711 integer :: ii, jj, nKS, iKS, nK, nS, iK, iS, fdUnit 1712 type(lnParams) :: params 1713 1714 call get_params(negf, params) 1715 1716 nKS = size(groupKS, dim=2) 1717 nK = size(kpoints, dim=2) 1718 nS = nKS/nK 1719 1720 do jj=1,size(matTot, dim=2) 1721 open(newunit=fdUnit, file=trim(regionLabels(jj))//'.dat') 1722 write(fdUnit,"(A)")'# Energy / eV States / e' 1723 do ii=1,size(matTot, dim=1) 1724 write(fdUnit,'(F12.6)',ADVANCE='NO') (params%Emin+(ii-1)*params%Estep) * Hartree__eV 1725 write(fdUnit,'(ES20.8)') matTot(ii,jj) 1726 enddo 1727 close(fdUnit) 1728 enddo 1729 1730 if (allocated(matSKRes)) then 1731 if (nKS > 1) then 1732 do jj = 1, size(matSKRes(:,:,:), dim=2) 1733 open(newunit=fdUnit, file=trim(regionLabels(jj))//'_kpoints.dat') 1734 write(fdUnit,"(A,I0)")'# NKpoints = ', nK 1735 write(fdUnit,"(A,I1)")'# NSpin = ', nS 1736 write(fdUnit,"(A)")'# <spin k1 k2 k3 weight> ' 1737 write(fdUnit,'(A1)', ADVANCE='NO') '# ' 1738 do iKS = 1, nKS 1739 iK = groupKS(1,iKS) 1740 iS = groupKS(2,iKS) 1741 write(fdUnit,'(i5.1)', ADVANCE='NO') iS 1742 write(fdUnit,'(es15.5, es15.5, es15.5, es15.5)', ADVANCE='NO') kpoints(:,iK),& 1743 & kWeights(iK) 1744 end do 1745 write(fdUnit,*) 1746 1747 do ii = 1, size(matSKRes(:,:,:), dim=1) 1748 write(fdUnit,'(f20.6)',ADVANCE='NO') (params%Emin+(ii-1)*params%Estep) * Hartree__eV 1749 do iKS = 1,nKS 1750 write(fdUnit,'(es20.8)',ADVANCE='NO') matSKRes(ii,jj, iKS) 1751 enddo 1752 write(fdUnit,*) 1753 enddo 1754 close(fdUnit) 1755 enddo 1756 end if 1757 end if 1758 1759 end subroutine write_files 1760 1761 1762 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1763 ! THIS is a first version of local current computation. 1764 ! It has been placed here since it depends on internal representations of DFTB 1765 ! 1766 ! NOTE: Limited to non-periodic systems s !!!!!!!!!!! 1767 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1768#:if WITH_MPI 1769 subroutine local_currents(mpicomm, groupKS, ham, over, & 1770 & neighbourList, nNeighbour, skCutoff, iAtomStart, iPair, img2CentCell, iCellVec, & 1771 & cellVec, rCellVec, orb, kPoints, kWeights, coord0, species0, speciesName, chempot, & 1772 & testArray) 1773 1774 !> MPI communicator 1775 type(mpifx_comm), intent(in) :: mpicomm 1776#:else 1777 subroutine local_currents(groupKS, ham, over, & 1778 & neighbourList, nNeighbour, skCutoff, iAtomStart, iPair, img2CentCell, iCellVec, & 1779 & cellVec, rCellVec, orb, kPoints, kWeights, coord0, species0, speciesName, chempot, & 1780 & testArray) 1781#:endif 1782 1783 !> kpoint and spin descriptor 1784 integer, intent(in) :: groupKS(:,:) 1785 1786 !> Hamiltonian and Overlap matrices 1787 real(dp), intent(in) :: ham(:,:), over(:) 1788 1789 !> Neighbor list container 1790 type(TNeighbourList), intent(in) :: neighbourList 1791 1792 !> nNeighbor list for SK interactions 1793 integer, intent(in) :: nNeighbour(:) 1794 1795 !> SK interaction cutoff 1796 real(dp), intent(in) :: skCutoff 1797 1798 !> Indices of staring atom block and Pairs 1799 integer, intent(in) :: iAtomStart(:), iPair(0:,:) 1800 1801 !> map of atoms to central cell 1802 integer, intent(in) :: img2CentCell(:), iCellVec(:) 1803 1804 !> translation vectors to lattice cells in units of lattice constants 1805 real(dp), allocatable, intent(in) :: cellVec(:,:) 1806 1807 !> Vectors to unit cells in absolute units 1808 real(dp), allocatable, intent(in) :: rCellVec(:,:) 1809 1810 !> Orbital descriptor 1811 type(TOrbitals), intent(in) :: orb 1812 1813 !> k-points and weights 1814 real(dp), intent(in) :: kPoints(:,:), kWeights(:) 1815 1816 !> central cell coordinates (folded to central cell) 1817 real(dp), intent(in) :: coord0(:,:) 1818 1819 !> Species indices for atoms in central cell 1820 integer, intent(in) :: species0(:) 1821 1822 !> Species Names (as in gen file) 1823 character(*), intent(in) :: speciesName(:) 1824 1825 ! We need this now for different fermi levels in colinear spin 1826 ! Note: spin polarized does not work with 1827 ! built-in potential (the unpolarized does) in the poisson 1828 ! I do not set the fermi because it seems that in libnegf it is 1829 ! not really needed 1830 real(dp), intent(in) :: chempot(:,:) 1831 1832 !> Array passed back to main for autotests (will become the output) 1833 real(dp), allocatable, intent(out) :: testArray(:,:) 1834 1835 1836 ! Local stuff --------------------------------------------------------- 1837 integer :: n0, nn, mm, mu, nu, nAtom, irow, nrow, ncont 1838 integer :: nKS, nKPoint, nSpin, iK, iKS, iS, inn, startn, endn, morb 1839 real(dp), dimension(:,:,:), allocatable :: lcurr 1840 real(dp) :: Im 1841 type(TNeighbourList) :: lc_neigh 1842 integer, dimension(:), allocatable :: lc_img2CentCell, lc_iCellVec, lc_species 1843 real(dp), dimension(:,:), allocatable :: lc_coord 1844 integer :: lc_nAllAtom 1845 real(dp) :: cutoff 1846 integer, parameter :: nInitNeigh=40 1847 complex(dp) :: c1,c2 1848 character(3) :: skp 1849 integer :: iSCCiter 1850 type(z_CSR), target :: csrDens, csrEDens 1851 type(z_CSR), pointer :: pCsrDens, pCsrEDens 1852 type(lnParams) :: params 1853 integer :: fdUnit 1854 1855 pCsrDens => csrDens 1856 pCsrEDens => csrEDens 1857 1858#:if WITH_MPI 1859 call negf_mpi_init(mpicomm) 1860#:endif 1861 call get_params(negf, params) 1862 1863 !Decide what to do with surface GFs. 1864 !sets readOldSGF: if it is 0 or 1 it is left so 1865 if (negf%readOldDM_SGFs.eq.COMPSAVE_SGF) then 1866 if(iSCCIter.eq.1) then 1867 call set_readOldDMsgf(negf, COMPSAVE_SGF) ! compute and save SGF on files 1868 else 1869 call set_readOldDMsgf(negf, READ_SGF) ! read from files 1870 endif 1871 endif 1872 1873 write(stdOut, *) 1874 write(stdOut, '(80("="))') 1875 write(stdOut, *) ' COMPUTING LOCAL CURRENTS ' 1876 write(stdOut, '(80("="))') 1877 write(stdOut, *) 1878 1879 nKS = size(groupKS, dim=2) 1880 nKPoint = size(kPoints, dim=2) 1881 nSpin = size(ham, dim=2) 1882 nAtom = size(orb%nOrbAtom) 1883 if (nKPoint > 999) then 1884 call error("Too many kpoints > 999 for local currents") 1885 end if 1886 1887 ! Create a symmetrized neighbour list extended to periodic cell in lc_coord 1888 if (any(iCellVec.ne.1)) then 1889 lc_nAllAtom = int((real(nAtom, dp)**(1.0_dp/3.0_dp) + 3.0_dp)**3) 1890 else 1891 lc_nAllAtom = nAtom 1892 end if 1893 allocate(lc_coord(3, lc_nAllAtom)) 1894 allocate(lc_species(lc_nAllAtom)) 1895 allocate(lc_img2CentCell(lc_nAllAtom)) 1896 allocate(lc_iCellVec(lc_nAllAtom)) 1897 call init(lc_neigh, nAtom, nInitNeigh) 1898 1899 call updateNeighbourListAndSpecies(lc_coord, lc_species, lc_img2CentCell, lc_iCellVec, & 1900 & lc_neigh, lc_nAllAtom, coord0, species0, skCutoff, rCellVec, symmetric=.true.) 1901 1902 allocate(lcurr(maxval(lc_neigh%nNeighbour), nAtom, nSpin)) 1903 lcurr = 0.0_dp 1904 1905 ! loop on k-points and spin 1906 do iKS = 1, nKS 1907 iK = groupKS(1, iKS) 1908 iS = groupKS(2, iKS) 1909 1910 write(stdOut,*) 'k-point',iK,'Spin',iS 1911 1912 ! We need to recompute Rho and RhoE ..... 1913 call foldToCSR(csrHam, ham(:,iS), kPoints(:,iK), iAtomStart, iPair, neighbourList%iNeighbour,& 1914 & nNeighbour, img2CentCell, iCellVec, CellVec, orb) 1915 call foldToCSR(csrOver, over, kPoints(:,iK), iAtomStart, iPair, neighbourList%iNeighbour, & 1916 & nNeighbour, img2CentCell, iCellVec, CellVec, orb) 1917 1918 call negf_density(iSCCIter, iS, iKS, pCsrHam, pCsrOver, chempot(:,iS), DensMat=pCsrDens) 1919 1920 ! Unless SGFs are not stored, read them from file 1921 if (negf%readOldDM_SGFs.ne.COMP_SGF) then 1922 call set_readOldDMsgf(negf, READ_SGF) 1923 end if 1924 1925 call negf_density(iSCCIter, iS, iKS, pCsrHam, pCsrOver, chempot(:,iS), EnMat=pCsrEDens) 1926 1927#:if WITH_MPI 1928 call mpifx_allreduceip(mpicomm, csrDens%nzval, MPI_SUM) 1929 call mpifx_allreduceip(mpicomm, csrEDens%nzval, MPI_SUM) 1930#:endif 1931 1932 write(skp,'(I3.3)') iK 1933 open(newUnit = fdUnit, file = 'lcurrents_'//skp//"_"//spin2ch(iS)//'.dat') 1934 1935 ! loop on central cell atoms and write local currents to all other 1936 ! interacting atoms within the cell and neighbour cells 1937 do mm = 1, nAtom 1938 1939 mOrb = orb%nOrbAtom(mm) 1940 iRow = iAtomStart(mm) 1941 1942 write(fdUnit,'(I5,3(F12.6),I4)',advance='NO') mm, lc_coord(:,mm), lc_neigh%nNeighbour(mm) 1943 1944 do inn = 1, lc_neigh%nNeighbour(mm) 1945 nn = lc_neigh%iNeighbour(inn, mm) 1946 n0 = lc_img2CentCell(nn) 1947 startn = iAtomStart(n0) 1948 endn = startn + orb%nOrbAtom(n0) - 1 1949 Im = 0.0_dp 1950 ! tracing orbitals of atoms n m 1951 ! More efficient without getel ? 1952 do mu = iRow, iRow+mOrb-1 1953 do nu = startn, endn 1954 c1 = conjg(getel(csrDens,mu,nu)) 1955 c2 = conjg(getel(csrEDens,mu,nu)) 1956 Im = Im + aimag(getel(csrHam,mu,nu)*c1 - getel(csrOver,mu,nu)*c2) 1957 enddo 1958 enddo 1959 ! pi-factor comes from Gn = rho * pi 1960 Im = Im * 2.0_dp*params%g_spin*pi*eovh*kWeights(iK) 1961 write(fdUnit,'(I5,ES17.8)',advance='NO') nn, Im 1962 lcurr(inn, mm, iS) = lcurr(inn, mm, iS) + Im 1963 enddo 1964 1965 write(fdUnit,*) 1966 enddo 1967 1968 close(fdUnit) 1969 1970 call destruct(csrDens) 1971 call destruct(csrEDens) 1972 1973 enddo 1974 1975 allocate(testArray(maxval(lc_neigh%nNeighbour),nAtom*nSpin)) 1976 testArray=0.0_dp 1977 ! Write the total current per spin channel 1978 do iS = 1, nSpin 1979 open(newUnit = fdUnit, file = 'lcurrents_'//spin2ch(iS)//'.dat') 1980 do mm = 1, nAtom 1981 write(fdUnit,'(I5,3(F12.6),I4)',advance='NO') mm, lc_coord(:,mm), lc_neigh%nNeighbour(mm) 1982 do inn = 1, lc_neigh%nNeighbour(mm) 1983 write(fdUnit,'(I5,ES17.8)',advance='NO') lc_neigh%iNeighbour(inn, mm), lcurr(inn,mm,iS) 1984 testArray(inn,(iS-1)*nAtom+mm) = lcurr(inn,mm,iS) 1985 end do 1986 write(fdUnit,*) 1987 end do 1988 end do 1989 close(fdUnit) 1990 deallocate(lcurr) 1991 1992 write(stdOut,*) 1993 call writeXYZFormat("supercell.xyz", lc_coord, lc_species, speciesName) 1994 write(stdOut,*) " <<< supercell.xyz written on file" 1995 1996 contains 1997 1998 !> labels from spin channel number 1999 pure function spin2ch(iS) result(ch) 2000 2001 !> channel number 2002 integer, intent(in) :: iS 2003 2004 !> label 2005 character(1) :: ch 2006 2007 character(1), parameter :: labels(2) = ['u', 'd'] 2008 2009 ch = labels(iS) 2010 2011 end function spin2ch 2012 2013 end subroutine local_currents 2014 2015 2016 !> pack dense matrices into CSR format 2017 subroutine MakeHHSS(H_all, S_all, HH, SS) 2018 2019 !> hamitonian matrix 2020 real(dp), intent(in) :: H_all(:,:) 2021 2022 !> overlap matrix 2023 real(dp), intent(in) :: S_all(:,:) 2024 2025 !> hamiltonian in CSR 2026 type(z_CSR), intent(inout) :: HH 2027 2028 !> overlap in CSR 2029 type(z_CSR), intent(inout) :: SS 2030 2031 integer :: i,j,k,l,m,n,NumStates, nnz 2032 2033 NumStates = negf%NumStates 2034 2035 nnz=0 2036 do i=1,NumStates 2037 do j=1,NumStates 2038 if ((i.eq.j).or.(abs(H_all(i,j)).gt.0.00001_dp)) then 2039 nnz = nnz+1 2040 end if 2041 end do 2042 end do 2043 2044 call destroy(HH) 2045 call create(HH, NumStates, NumStates, nnz) 2046 call destroy(SS) 2047 call create(SS, NumStates, NumStates, nnz) 2048 2049 HH%rowpnt(1)=1 2050 nnz=0 2051 do i=1,NumStates 2052 k=0 2053 do j=1,NumStates 2054 if((i.eq.j).or.(abs(H_all(i,j)).gt.0.00001_dp)) then 2055 k=k+1 2056 nnz=nnz+1 2057 if(i.eq.j) then 2058 HH%nzval(nnz)= H_all(i,j) 2059 SS%nzval(nnz)= S_all(i,j) 2060 else 2061 HH%nzval(nnz)= H_all(i,j) 2062 SS%nzval(nnz)= S_all(i,j) 2063 end if 2064 HH%colind(nnz)=j 2065 !print *, HH%nnz, negf%H_all(i,j), HH%nzval(HH%nnz), SS%nzval(HH%nnz) !debug 2066 end if 2067 end do 2068 HH%rowpnt(i+1)=HH%rowpnt(i)+k 2069 end do 2070 2071 SS%colind=HH%colind 2072 SS%rowpnt=HH%rowpnt 2073 2074 end subroutine MakeHHSS 2075 2076 2077 !> form orthogonal matrices via Lowdin transform for whole system 2078 subroutine orthogonalization(H,S) 2079 2080 !> hamitonian matrix 2081 real(dp), intent(inout) :: H(:,:) 2082 2083 !> overlap matrix 2084 real(dp), intent(inout) :: S(:,:) 2085 2086 integer :: i,m,n1_first,n1_last,n2_first,n2_last 2087 integer :: INFO, N 2088 real(dp), allocatable :: A(:,:), WORK(:), W(:) 2089 real(dp), allocatable :: B(:,:),C(:,:) 2090 2091 N = size(H, dim=1) 2092 2093 if (N /= negf%NumStates) then 2094 call error('orthogonalization: negf init NumStates error') 2095 end if 2096 2097 allocate(A(N,N),WORK(3*N),W(N)) 2098 allocate(B(N,N),C(N,N)) 2099 W=0.0_dp 2100 WORK=0.0_dp 2101 A=0.0_dp 2102 B=0.0_dp 2103 C=0.0_dp 2104 2105 A=S 2106 2107 call DSYEV('V','U',N,A,N,W,WORK,3*N,INFO ) 2108 2109 !print *,'U matrix, Eigenvectors for S diagonalization' 2110 !do i=1,N 2111 ! write(*,*)A(i,1:N) 2112 !end do 2113 2114 !print *,'U matrix unitarity check' 2115 !B=matmul(transpose(A),A) 2116 !do i=1,N 2117 ! write(*,*)B(i,1:N) 2118 !end do 2119 2120 B(:,:) = matmul(transpose(A), matmul(S, A)) 2121 2122 do i=1,N 2123 B(i,i) = 1.0_dp / sqrt(B(i,i)) 2124 end do 2125 2126 !C=sqrt(S) 2127 C(:,:) = matmul(A, matmul(B, transpose(A))) 2128 2129 !print *,'sqrt(S) inverted' 2130 !do i=1,N 2131 ! write(*,*) C(i,1:N) 2132 !end do 2133 2134 !print *,'S unity check' 2135 !B=matmul(transpose(C),matmul(S,C)) 2136 !do i=1,N 2137 ! write(*,*) B(i,1:N) 2138 !end do 2139 2140 !print *,'H_dftb before orthogonalization' 2141 !do i=1,N 2142 ! write(*,*) H(i,1:N) 2143 !end do 2144 2145 H(:,:) = matmul(transpose(C), matmul(H, C)) 2146 2147 !print *,'H_dftb_orth before replacement' 2148 !do i=1,N 2149 ! write(*,*) H(i,1:N) 2150 !end do 2151 2152 ! COPY THE FIRST CONTACT PL ONTO THE SECOND 2153 do m = 1, negf%str%num_conts 2154 n1_first = negf%str%mat_B_start(m) 2155 n1_last = (negf%str%mat_C_end(m)+negf%str%mat_B_start(m))/2 2156 n2_first = n1_last + 1 2157 n2_last = negf%str%mat_C_end(m) 2158 !print *,n1_first,n1_last,n2_first,n2_last 2159 H(n2_first:n2_last,n2_first:n2_last) = H(n1_first:n1_last,n1_first:n1_last) 2160 end do 2161 2162 !print *,'H_dftb_orth after replacement' 2163 !do i=1,N 2164 ! write(*,*) H(i,1:N) 2165 !end do 2166 2167 S = 0.0_dp 2168 do i=1,N 2169 S(i,i) = 1.0_dp 2170 end do 2171 2172 !Save H_dftb_orth.mtr to file 2173 !open(12,file='H_dftb_orth.mtr',action="write") 2174 !do i = 1,N 2175 ! write(12,*) H(i,1:N)* Hartree__eV 2176 !end do 2177 !close(12) 2178 2179 end subroutine orthogonalization 2180 2181 2182 !> Orthogonalise basis in device region 2183 subroutine orthogonalization_dev(H, S) 2184 2185 !> hamilonian matrix 2186 real(dp), intent(inout) :: H(:,:) 2187 2188 !> overlap matrix 2189 real(dp), intent(inout) :: S(:,:) 2190 2191 2192 integer :: i,m,n1_first,n1_last,n2_first,n2_last 2193 integer :: INFO, N, N2 2194 real(dp), allocatable :: A(:,:), WORK(:), W(:) 2195 real(dp), allocatable :: B(:,:), U(:,:), C(:,:) 2196 2197 2198 N = size(H, dim=1) 2199 if (N /= negf%NumStates) then 2200 call error('orthogonalization: negf init NumStates error') 2201 end if 2202 2203 N2=negf%str%central_dim 2204 2205 allocate(A(N2,N2),WORK(3*N2),W(N2)) 2206 allocate(B(N2,N2), U(N2,N2)) 2207 W=0.0_dp 2208 WORK=0.0_dp 2209 2210 A = S(1:N2,1:N2) 2211 U = A 2212 2213 call DSYEV('V','U',N2,U,N2,W,WORK,3*N2,INFO ) 2214 2215 !print *,'U matrix, Eigenvectors for S diagonalization' 2216 !do i=1,N2 2217 ! write(*,*) U(i,1:N2) 2218 !end do 2219 2220 !U matrix unitarity check 2221 !B = matmul(transpose(U),U) 2222 !do i=1,N2 2223 ! if (abs(B(i,i)-1.0_dp) > 1.d-9 .or. any(abs(B(i,1:i-1)) > 1.d-9 ) then 2224 ! print*, 'ERROR: U is not unitary', B(i,:) 2225 ! stop 2226 ! end if 2227 !end do 2228 2229 B = matmul(transpose(U),matmul(A,U)) 2230 2231 do i=1,N2 2232 B(i,i)=1.0_dp/sqrt(B(i,i)) 2233 end do 2234 2235 2236 ! Now A = S^-1/2 2237 A = matmul(U,matmul(B,transpose(U))) 2238 !print *,'sqrt(S) inverted' 2239 !do i=1,N2 2240 ! write(*,*) A(i,1:N2) 2241 !end do 2242 2243 deallocate(U, B) 2244 2245 allocate(C(N,N)) 2246 C=0.0_dp 2247 do i = N2+1,N 2248 C(i,i)=1.0_dp 2249 end do 2250 C(1:N2,1:N2) = A 2251 2252 deallocate(A) 2253 2254 !C=sqrt(S) big matrix 2255 !print *,'C=sqrt(S) big matrix' 2256 !do i=1,N 2257 ! write(*,*) C(i,1:N) 2258 !end do 2259 2260 !print *,'H_dftb before orthogonalization' 2261 !do i=1,N 2262 ! write(*,*) H(i,1:N) 2263 !end do 2264 2265 H = matmul(transpose(C),matmul(H,C)) 2266 S = matmul(transpose(C),matmul(S,C)) 2267 2268 !print *,'H_dftb_orth before replacement' 2269 !do i=1,N 2270 ! write(*,*) H(i,1:N) 2271 !end do 2272 2273 ! COPY THE FIRST CONTACT PL ONTO THE SECOND 2274 do m=1,negf%str%num_conts 2275 n1_first=negf%str%mat_B_start(m) 2276 n1_last=(negf%str%mat_C_end(m)+negf%str%mat_B_start(m))/2 2277 n2_first=n1_last+1 2278 n2_last=negf%str%mat_C_end(m) 2279 !print *,n1_first,n1_last,n2_first,n2_last 2280 H(n2_first:n2_last,n2_first:n2_last) = H(n1_first:n1_last,n1_first:n1_last) 2281 S(n2_first:n2_last,n2_first:n2_last) = S(n1_first:n1_last,n1_first:n1_last) 2282 end do 2283 2284 !print *,'H_dftb_orth after replacement' 2285 !do i=1,N 2286 ! write(*,*) H(i,1:N) 2287 !end do 2288 2289 !print *,'S_dftb_orth after replacement' 2290 !do i=1,N 2291 ! write(*,*) S(i,1:N) 2292 !end do 2293 2294 !Save H_dftb_orth.mtr to file 2295 !open(12,file='H_dftb_orth.mtr',action="write") 2296 !do i=1,N 2297 ! write(12,*) H(i,1:N)*Hartree__eV 2298 !end do 2299 !close(12) 2300 2301 !Save S_dftb_orth.mtr to file 2302 !open(12,file='S_dftb_orth.mtr',action="write") 2303 !do i=1,N 2304 ! write(12,*) S(i,1:N) 2305 !end do 2306 !close(12) 2307 2308 end subroutine orthogonalization_dev 2309 2310 2311end module negf_int 2312