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, 2012, nickpapior@gmail.com 10! 11module m_hs_matrix 12! 13! Routines that are used for the distribution of the Hamiltonian 14! and scattering matrices into a full size matrix instead of sparse matrices. 15! It has several options for creating different kinds of matrices. 16! This module is a serial version. It requires that the Node has the full 17! sparse matrix available. 18! 19! Specifically it can be used to remove the z-direction connections. 20! Also the inner-cell distances are in the xij arrays. Through routine calls 21! these inner cell distances can be removed. 22! 23! The reason for having this is future use of routines which can utilize the 24! Hamiltonian created in different ways. 25! Say you want to test calculate a transmission from a Hamiltonian which is 26! created by SIESTA. For that purpose you need to remove the cell connection 27! in the z-direction. 28! 29! The usage of this module is highly encouraged in future utilities where 30! the need for the Hamiltonian and/or overlap matrix is needed. 31! With this in mind the code can further be optimized for speed. 32! However, for normal sizes it is quite fast. 33! 34! Also for testing against Inelastica the use of inner-cell distances is not 35! used. Therefore the Hamiltonians can not be numerically compared. 36! If this is to be enforced later on. The option is there. 37! 38! The use of this module is a straight forward call: 39! 40! call set_HS_matrix(Gamma,ucell,na_u,no_u,no_s,maxnh, & 41! xij,numh,listhptr,listh,H,S, & 42! k,Hk,Sk, & 43! xa,iaorb, & 44! RemZConnection,RemUCellDistances,RemNFirstOrbitals,RemNLastOrbitals) 45! 46! xa, iaorb, RemZConnection, RemUCellDistances, RemNFirstOrbitals, RemNLastOrbitals 47! are all optional arguments. xa and iaorb are required if RemZConnection or RemUCellDistances 48! are true. 49! Notice, that any calls to the optional arguments MUST be with keywords! Otherwise 50! the program will end! 51! We have added so that when iaorb is required, you could also use lasto. This 52! makes it more intuitive as lasto is often more accesible. 53! It on the other hand has a small overhead when using lasto (negligeble). 54! 55! Gamma denotes whether it is a Gamma calculation (if true, it will not 56! add k-phases, no matter if k /= \Gamma-point. 57! na_u,no_u,no_s,maxnh,xij,numh,listhptr,listh,H,S are all variables 58! needed in the definition of the entire H and S matrices in the sparse format. 59! 60! k is the k-point that will be created for the Hamiltonian. 61! Hk and Sk are returned for the user. 62! 63! The RemZConnection can be used to remove any matrix elements <i|H|j> where i and j 64! are connections in the next unit cell in the Z-direction. 65! This is usefull if one wishes to see the matrix as it would look while doing 66! transmission calculations. 67! 68! RemUCellDistances can be set to true so that inner-cell distances are removed. 69! Several people on the SIESTA mailing list have "complained/asked questions" about 70! the difference in the Hamiltonians which are not always constructed similarly. 71! However, inner cell differences can be neglected. This is merely to get the same 72! matrix as is created in Inelastica, for example. 73! 74! RemNFirstOrbitals can be used to fully remove that many states from the start of the 75! Hamiltonian. This is used when there are buffer regions which should be disregarded. 76! RemNLastOrbitals is the same, albeit in the end of the Hamiltonian. 77! 78! If the sizes of the incoming pointers Hk and Sk do not match the above 79! they will be reallocated to the correct size (without notifying the user). 80! 81! Also there is a routine for obtaining an arbitrary transfer matrix. 82! The interface is very much the same as that for the set_HS_matrix. 83! However, here it has this interface: 84! 85! call set_HS_transfermatrix(Gamma,ucell,na_u,no_u,no_s,maxnh, & 86! xij,numh,listhptr,listh,H,S, & 87! k,transfer_cell,HkT,SkT,xa,iaorb, & 88! RemUCellDistances,RemNFirstOrbitals,RemNLastOrbitals) 89! 90! 'transfer_cell' is an integer vector of dimen 3 with each indices denoting the transfer matrix in that 91! direction. 92! For instance taking: 93! transfer_cell(:) = (/ 1 , 3 , 1 /) 94! Will generate the transfer cell to the first neighbouring 'x', third neighbouring 'y' and first 95! neighbouring 'z' cell. 96! Often it can be necessary to have access to only the transfer matrix in one direction. 97! In such cases you can supply: 98! transfer_cell(:) = (/ TRANSFER_ALL , TRANSFER_ALL , 2 /) 99! The TRANSFER_ALL is a variable in this module and tells that it should 100! use all possible transfer cells in that direction. 101! 102! In this conjunction the routine 103! call set_HS_available_transfers(Gamma,ucell,na_u,no_u,no_s,maxnh, & 104! xij,numh,listhptr,listh,xa,iaorb,transfer_cell) 105! is invaluable. It returns in transfer_cell the allowed transfer cell units for 106! the system. Notice that transfer_cell is a 2x3 matrix with the formatting like this: 107! 108! transfer_cell(1,1) : minumum transfer cell in the x-direction 109! transfer_cell(2,1) : maximum transfer cell in the x-direction 110! transfer_cell(1,2) : minumum transfer cell in the y-direction 111! transfer_cell(2,2) : maximum transfer cell in the y-direction 112! transfer_cell(1,3) : minumum transfer cell in the z-direction 113! transfer_cell(2,3) : maximum transfer cell in the z-direction 114! 115! Usually the following holds (for obvious reasons): 116! transfer_cell(1,i) = -transfer_cell(2,i) 117! This routine can thus be used to retrieve the allowed transfer cells before using 118! set_HS_transfermatrix. 119! 120! Notice that RemZConnection is NOT available (it doesn't make sense). 121! Furthermore, xa and iaorb are needed in order to determine the transfer cell. 122! Thus they are no longer optional but mandatory arguments. 123! 124! 125! Furthermore there are the routines: 126! call matrix_rem_left_right(no_tot,Hk,Sk,no_L,no_R) 127! and 128! call matrix_symmetrize(no_tot,Hk,Sk,Ef) 129! Which are used to remove connections of left/right regions 130! and used to symmetrize and shift the Hamiltonian Ef, respectively. 131! 132! NOTICE that a call to matrix_symmetrize is almost always needed! 133! EVEN in the case of the transfer matrix. 134 135 use precision, only : dp 136 137 implicit none 138 139 private 140 141 interface set_HS_matrix 142 module procedure set_HS_matrix_1d 143 module procedure set_HS_matrix_2d 144 end interface set_HS_matrix 145 146 interface set_HS_transfermatrix 147 module procedure set_HS_transfermatrix_1d 148 module procedure set_HS_transfermatrix_2d 149 end interface set_HS_transfermatrix 150 151 interface matrix_rem_left_right 152 module procedure matrix_rem_left_right_1d 153 module procedure matrix_rem_left_right_2d 154 end interface matrix_rem_left_right 155 156 interface matrix_symmetrize 157 module procedure matrix_symmetrize_1d 158 module procedure matrix_symmetrize_2d 159 end interface matrix_symmetrize 160 161 public :: set_HS_matrix 162 public :: set_HS_transfermatrix 163 public :: matrix_rem_left_right 164 public :: matrix_symmetrize 165 public :: set_HS_available_transfers 166 167 integer, parameter, public :: TRANSFER_ALL = -999999 168 169contains 170 171!***************** 172! Setting the Hamiltonian for a specific k-point. 173! It requires that the Node has the full sparse matrix available 174!***************** 175 subroutine set_HS_matrix_1d(Gamma,ucell,na_u,no_u,no_s,maxnh, & 176 xij,numh,listhptr,listh,H,S, & 177 k,Hk,Sk, & 178 DUMMY, & ! Ensures that the programmer makes EXPLICIT keywork passing 179 xa,iaorb,lasto, & 180 RemZConnection,RemUCellDistances,RemNFirstOrbitals,RemNLastOrbitals) 181 use sys, only : die 182 use alloc, only : re_alloc 183 use geom_helper, only : ucorb 184 use cellSubs, only : reclat 185 186! *********************** 187! * INPUT variables * 188! *********************** 189 logical, intent(in) :: Gamma ! Is it a Gamma Calculation? 190 real(dp), intent(in) :: ucell(3,3) ! The unit cell of system 191 integer, intent(in) :: na_u ! Unit cell atoms 192 integer, intent(in) :: no_u ! Unit cell orbitals 193 integer, intent(in) :: no_s ! Supercell orbitals 194 integer, intent(in) :: maxnh ! Hamiltonian size 195 real(dp), intent(in) :: xij(3,maxnh) ! differences with unitcell, differences with unitcell 196 integer, intent(in) :: numh(no_u),listhptr(no_u) 197 integer, intent(in) :: listh(maxnh) 198 real(dp), intent(in) :: H(maxnh) ! Hamiltonian 199 real(dp), intent(in) :: S(maxnh) ! Overlap 200 real(dp), intent(in) :: k(3) ! k-point in [1/Bohr] 201! *********************** 202! * OUTPUT variables * 203! *********************** 204 complex(dp), pointer :: Hk(:), Sk(:) 205 206! *********************** 207! * OPTIONAL variables * 208! *********************** 209 logical, intent(in), optional :: DUMMY ! Do not supply this, it merely requires the coder 210! ! to use the keyworded arguments! 211 real(dp), intent(in),optional :: xa(3,na_u) ! Atomic coordinates (needed for RemZConnection & RemUCellDistances) 212 integer, intent(in), optional :: iaorb(no_u) ! The equivalent atomic index for a given orbital (needed for RemUCellDistances) 213 integer, intent(in), optional :: lasto(0:na_u) ! The number of orbitals on each atom (needed for RemUCellDistances) 214 logical, intent(in), optional :: RemZConnection, RemUCellDistances 215 integer, intent(in), optional :: RemNFirstOrbitals, RemNLastOrbitals 216 217! *********************** 218! * LOCAL variables * 219! *********************** 220 real(dp) :: recell(3,3) ! The reciprocal unit cell 221 real(dp) :: xo(3), xc 222 real(dp), allocatable :: xuo(:) 223 real(dp) :: kxij 224 complex(dp) :: cphase 225 integer :: no_tot 226 integer, allocatable :: liaorb(:) 227 integer :: i,j,iu,iuo,juo,iind,ind 228 logical :: l_RemZConnection, l_RemUCellDistances 229 integer :: l_RemNFirstOrbitals, l_RemNLastOrbitals 230 231 if ( present(DUMMY) ) & 232 call die("You must specify the keyworded arguments & 233 &for set_HS_matrix") 234 235 ! Option collecting 236 l_RemZConnection = .false. 237 if ( present(RemZConnection) ) & 238 l_RemZConnection = RemZConnection 239 if (l_RemZConnection .and. .not. present(xa)) & 240 call die("You need xa in set_HS_matrix when removing & 241 &the z-connection.") 242 if (l_RemZConnection .and. & 243 (.not. present(iaorb) .and. .not. present(lasto) ) ) & 244 call die("You need iaorb or lasto in set_HS_matrix when removing & 245 &the z-connection.") 246 l_RemUCellDistances = .false. 247 if ( present(RemUCellDistances) ) & 248 l_RemUCellDistances = RemUCellDistances 249 if (l_RemUCellDistances .and. .not. present(xa)) & 250 call die("You need xa in set_HS_matrix when removing & 251 &unit cell distances.") 252 if (l_RemUCellDistances .and. & 253 (.not. present(iaorb) .and. .not. present(lasto) ) ) & 254 call die("You need iaorb or lasto in set_HS_matrix when removing & 255 &unit cell distances.") 256 257 ! Make l_RemNFirstOrbitals contain the number of orbitals 258 ! to be removed from the Hamiltonian in the BEGINNING 259 l_RemNFirstOrbitals = 0 260 if ( present(RemNFirstOrbitals) ) & 261 l_RemNFirstOrbitals = RemNFirstOrbitals 262 ! Make l_RemNLastOrbitals contain the number of orbitals 263 ! to be removed from the Hamiltonian in the END 264 l_RemNLastOrbitals = 0 265 if ( present(RemNLastOrbitals) ) & 266 l_RemNLastOrbitals = RemNLastOrbitals 267 no_tot = no_u - (l_RemNLastOrbitals + l_RemNFirstOrbitals) 268 269 270 call re_alloc(Hk,1,no_tot*no_tot,name='Hk',routine='set_HS') 271 call re_alloc(Sk,1,no_tot*no_tot,name='Sk',routine='set_HS') 272 273 if ( l_RemZConnection ) then 274 ! Prepare the cell to calculate the index of the atom 275 call reclat(ucell,recell,0) ! Without 2*Pi 276 277 ! Find the actual coordinates of the orbitals in the form of the sparse matrices 278 ! Notice that this array is without the removed orbitals 279 allocate(xuo(no_tot)) 280 call memory('A','D',no_tot,'set_HS') 281 282 if ( present(iaorb) ) then 283 do iuo = 1 , no_tot 284 i = iaorb(iuo + l_RemNFirstOrbitals) 285 xuo(iuo) = & 286 xa(1,i) * recell(1,3) + & 287 xa(2,i) * recell(2,3) + & 288 xa(3,i) * recell(3,3) 289 end do !io in uc 290 else if ( present(lasto) ) then 291 do j = 1 , na_u 292 do i = lasto(j-1) + 1 , lasto(j) 293 if ( i <= l_RemNFirstOrbitals ) cycle 294 if ( no_tot < i - l_RemNFirstOrbitals ) cycle 295 xuo(i-l_RemNFirstOrbitals) = & 296 xa(1,j) * recell(1,3) + & 297 xa(2,j) * recell(2,3) + & 298 xa(3,j) * recell(3,3) 299 end do !io in uc 300 end do 301 end if 302 303 end if 304 305 ! Create the orb => atom index array 306 if ( l_RemUCellDistances ) then 307 allocate(liaorb(no_tot)) 308 call memory('A','I',no_tot,'set_HS') 309 310 if ( present(iaorb) ) then 311 do iuo = 1 , no_tot 312 liaorb(iuo) = iaorb(iuo + l_RemNFirstOrbitals) 313 end do !io in uc 314 else if ( present(lasto) ) then 315 ind = 0 316 do j = 1 , na_u 317 do i = lasto(j-1) + 1 , lasto(j) 318 if ( i <= l_RemNFirstOrbitals ) cycle 319 if ( no_tot < i - l_RemNFirstOrbitals ) cycle 320 ind = ind + 1 321 liaorb(ind) = j 322 end do !io in uc 323 end do 324 end if 325 326 end if 327 328! 329! Setup H,S for this k-point: 330! 331 do i = 1,no_tot*no_tot 332 Hk(i) = dcmplx(0.d0,0.d0) 333 Sk(i) = dcmplx(0.d0,0.d0) 334 end do 335 336 xo(:) = 0.0_dp 337 338 setup_HS: if (.not.Gamma ) then 339 340 do iuo = 1 , no_tot 341 iu = iuo + l_RemNFirstOrbitals 342 iind = listhptr(iu) 343 do j = 1,numh(iu) 344 ind = iind + j 345 juo = ucorb(listh(ind),no_u) - l_RemNFirstOrbitals 346 347 ! Cycle if we are not in the middle region 348 if ( juo < 1 .or. no_tot < juo ) cycle 349 350 ! If we have a removal of the z-direction do the following 351 if ( l_RemZConnection ) then 352 xc = - (xuo(juo) - xuo(iuo)) 353 xc = xc + xij(1,ind) * recell(1,3) + & 354 xij(2,ind) * recell(2,3) + & 355 xij(3,ind) * recell(3,3) 356 if ( nint(xc) /= 0 ) cycle 357 end if 358 359 ! We also wish to remove the connection in 360 ! in the inner cell 361 ! I suspect we have a problem here 362 ! We may need a check for being in the unitcell 363 ! i.e.: 364 ! if ( l_RemUCellDistances ) then 365 ! if ( listh(ind) > no_u ) then ! As we check in the unit cell! 366 ! xo(1) = xa(1,liaorb(juo)) & 367 ! - xa(1,liaorb(iuo)) 368 ! xo(2) = xa(2,liaorb(juo)) & 369 ! - xa(2,liaorb(iuo)) 370 ! xo(3) = xa(3,liaorb(juo)) & 371 ! - xa(3,liaorb(iuo)) 372 ! else 373 ! xo = 0.0_dp 374 ! end if 375 ! end if 376 if ( l_RemUCellDistances ) then 377 xo(1) = xa(1,liaorb(juo)) & 378 - xa(1,liaorb(iuo)) 379 xo(2) = xa(2,liaorb(juo)) & 380 - xa(2,liaorb(iuo)) 381 xo(3) = xa(3,liaorb(juo)) & 382 - xa(3,liaorb(iuo)) 383 end if 384 385 kxij = & 386 k(1) * (xij(1,ind) - xo(1)) + & 387 k(2) * (xij(2,ind) - xo(2)) + & 388 k(3) * (xij(3,ind) - xo(3)) 389 cphase = exp(dcmplx(0d0,kxij)) 390 i = iuo+(juo-1)*no_tot 391 Hk(i) = Hk(i)+H(ind)*cphase 392 Sk(i) = Sk(i)+S(ind)*cphase 393 end do 394 end do 395 396 else setup_HS 397 ! It is not a Gamma calculation, thus we do not have any 398 ! neighbouring cells etc. 399 do iuo = 1 , no_tot 400 iu = iuo + l_RemNFirstOrbitals 401 iind = listhptr(iu) 402 do j = 1,numh(iu) 403 ind = iind + j 404 juo = listh(ind) - l_RemNFirstOrbitals 405 406 ! Cycle if we are not in the middle region 407 if ( juo < 1 .or. no_tot < juo ) cycle 408 409 ! If we have a removal of the z-direction do the following 410 if ( l_RemZConnection ) then 411 xc = - (xuo(juo) - xuo(iuo)) 412 xc = xc + xij(1,ind) * recell(1,3) + & 413 xij(2,ind) * recell(2,3) + & 414 xij(3,ind) * recell(3,3) 415 if ( nint(xc) /= 0 ) cycle 416 end if 417 418 i = iuo+(juo-1)*no_tot 419 Hk(i) = Hk(i)+H(ind) 420 Sk(i) = Sk(i)+S(ind) 421 end do 422 end do 423 424 end if setup_HS 425 426 if ( l_RemUCellDistances ) then 427 call memory('D','I',no_tot,'set_HS') 428 deallocate(liaorb) 429 end if 430 431 if ( l_RemZConnection ) then 432 call memory('D','D',no_tot,'set_HS') 433 deallocate(xuo) 434 end if 435 436 end subroutine set_HS_matrix_1d 437 438!***************** 439! Setting the Hamiltonian for a specific k-point. 440! It requires that the Node has the full sparse matrix available 441!***************** 442 subroutine set_HS_matrix_2d(Gamma,ucell,na_u,no_u,no_s,maxnh, & 443 xij,numh,listhptr,listh,H,S, & 444 k,Hk,Sk, & 445 DUMMY, & ! Ensures that the programmer makes EXPLICIT keywork passing 446 xa,iaorb,lasto, & 447 RemZConnection,RemUCellDistances,RemNFirstOrbitals,RemNLastOrbitals) 448 use sys, only : die 449 use alloc, only : re_alloc 450 use geom_helper, only : ucorb 451 use cellSubs, only : reclat 452 453! *********************** 454! * INPUT variables * 455! *********************** 456 logical, intent(in) :: Gamma ! Is it a Gamma Calculation? 457 real(dp), intent(in) :: ucell(3,3) ! The unit cell of system 458 integer, intent(in) :: na_u ! Unit cell atoms 459 integer, intent(in) :: no_u ! Unit cell orbitals 460 integer, intent(in) :: no_s ! Total orbitals 461 integer, intent(in) :: maxnh ! Hamiltonian size 462 real(dp), intent(in) :: xij(3,maxnh) ! differences with unitcell, differences with unitcell 463 integer, intent(in) :: numh(no_u),listhptr(no_u) 464 integer, intent(in) :: listh(maxnh) 465 real(dp), intent(in) :: H(maxnh) ! Hamiltonian 466 real(dp), intent(in) :: S(maxnh) ! Overlap 467 real(dp), intent(in) :: k(3) ! k-point in [1/Bohr] 468! *********************** 469! * OUTPUT variables * 470! *********************** 471 complex(dp), pointer :: Hk(:,:), Sk(:,:) 472! *********************** 473! * OPTIONAL variables * 474! *********************** 475 logical, intent(in), optional :: DUMMY ! Do not supply this, it merely requires the coder 476! ! to use the keyworded arguments! 477 real(dp), intent(in),optional :: xa(3,na_u) ! Atomic coordinates (needed for RemZConnection & RemUCellDistances) 478 integer, intent(in), optional :: iaorb(no_u) ! The equivalent atomic index for a given orbital (needed for RemUCellDistances) 479 integer, intent(in), optional :: lasto(0:na_u) ! The number of orbitals on each atom (needed for RemUCellDistances) 480 logical, intent(in), optional :: RemZConnection, RemUCellDistances 481 integer, intent(in), optional :: RemNFirstOrbitals, RemNLastOrbitals 482 483! *********************** 484! * LOCAL variables * 485! *********************** 486 real(dp) :: recell(3,3) 487 real(dp) :: xo(3), xc 488 real(dp), allocatable :: xuo(:) 489 integer :: no_tot 490 real(dp) :: kxij 491 complex(dp) :: cphase 492 integer, allocatable :: liaorb(:) 493 integer :: i,j,iuo,iu,juo,iind,ind 494 logical :: l_RemZConnection, l_RemUCellDistances 495 integer :: l_RemNFirstOrbitals, l_RemNLastOrbitals 496 497 if ( present(DUMMY) ) & 498 call die("You must specify the keyworded arguments & 499 &for set_HS_matrix") 500 501 ! Option collecting 502 l_RemZConnection = .false. 503 if ( present(RemZConnection) ) & 504 l_RemZConnection = RemZConnection 505 if (l_RemZConnection .and. .not. present(xa)) & 506 call die("You need xa in set_HS_matrix when removing & 507 &the z-connection.") 508 if (l_RemZConnection .and. & 509 (.not. present(iaorb) .and. .not. present(lasto) ) ) & 510 call die("You need iaorb or lasto in set_HS_matrix when removing & 511 &the z-connection.") 512 l_RemUCellDistances = .false. 513 if ( present(RemUCellDistances) ) & 514 l_RemUCellDistances = RemUCellDistances 515 if (l_RemUCellDistances .and. .not. present(xa)) & 516 call die("You need xa in set_HS_matrix when removing & 517 &unit cell distances.") 518 if (l_RemUCellDistances .and. & 519 (.not. present(iaorb) .and. .not. present(lasto) ) ) & 520 call die("You need iaorb or lasto in set_HS_matrix when removing & 521 &unit cell distances.") 522 523 524 ! Make l_RemNFirstOrbitals contain the number of orbitals 525 ! to be removed from the Hamiltonian in the BEGINNING 526 l_RemNFirstOrbitals = 0 527 if ( present(RemNFirstOrbitals) ) & 528 l_RemNFirstOrbitals = RemNFirstOrbitals 529 ! Make l_RemNLastOrbitals contain the number of orbitals 530 ! to be removed from the Hamiltonian in the END 531 l_RemNLastOrbitals = 0 532 if ( present(RemNLastOrbitals) ) & 533 l_RemNLastOrbitals = RemNLastOrbitals 534 no_tot = no_u - (l_RemNLastOrbitals + l_RemNFirstOrbitals) 535 536 537 call re_alloc(Hk,1,no_tot,1,no_tot,name='Hk',routine='set_HS') 538 call re_alloc(Sk,1,no_tot,1,no_tot,name='Sk',routine='set_HS') 539 540 if ( l_RemZConnection ) then 541 ! Prepare the cell to calculate the index of the atom 542 call reclat(ucell,recell,0) ! Without 2*Pi 543 544 ! Find the actual coordinates of the orbitals in the form of the sparse matrices 545 ! Notice that this array is without the removed orbitals 546 allocate(xuo(no_tot)) 547 call memory('A','D',no_tot,'set_HS') 548 549 if ( present(iaorb) ) then 550 do iuo = 1 , no_tot 551 i = iaorb(iuo + l_RemNFirstOrbitals) 552 xuo(iuo) = & 553 xa(1,i) * recell(1,3) + & 554 xa(2,i) * recell(2,3) + & 555 xa(3,i) * recell(3,3) 556 end do !io in uc 557 else if ( present(lasto) ) then 558 do j = 1 , na_u 559 do i = lasto(j-1) + 1 , lasto(j) 560 if ( i <= l_RemNFirstOrbitals ) cycle 561 if ( no_tot < i - l_RemNFirstOrbitals ) cycle 562 xuo(i-l_RemNFirstOrbitals) = & 563 xa(1,j) * recell(1,3) + & 564 xa(2,j) * recell(2,3) + & 565 xa(3,j) * recell(3,3) 566 end do !io in uc 567 end do 568 end if 569 570 end if 571 572 if ( l_RemUCellDistances ) then 573 allocate(liaorb(no_tot)) 574 call memory('A','I',no_tot,'set_HS') 575 576 if ( present(iaorb) ) then 577 do iuo = 1 , no_tot 578 liaorb(iuo) = iaorb(iuo + l_RemNFirstOrbitals) 579 end do !io in uc 580 else if ( present(lasto) ) then 581 ind = 0 582 do j = 1 , na_u 583 do i = lasto(j-1) + 1 , lasto(j) 584 if ( i <= l_RemNFirstOrbitals ) cycle 585 if ( no_tot < i - l_RemNFirstOrbitals ) cycle 586 ind = ind + 1 587 liaorb(ind) = j 588 end do !io in uc 589 end do 590 end if 591 592 end if 593 594 595! 596! Setup H,S for this k-point: 597! 598 do juo = 1,no_tot 599 do iuo = 1,no_tot 600 Hk(iuo,juo) = dcmplx(0.d0,0.d0) 601 Sk(iuo,juo) = dcmplx(0.d0,0.d0) 602 end do 603 end do 604 605 xo(:) = 0.0_dp 606 607 setup_HS: if (.not.Gamma ) then 608 609 do iuo = 1,no_tot 610 iu = iuo + l_RemNFirstOrbitals 611 iind = listhptr(iu) 612 do j = 1,numh(iu) 613 ind = iind + j 614 juo = ucorb(listh(ind),no_u) - l_RemNFirstOrbitals 615 616 ! Cycle if we are not in the middle region 617 if ( juo < 1 .or. no_tot < juo ) cycle 618 619 ! If we have a removal of the z-direction do the following 620 if ( l_RemZConnection ) then 621 xc = - (xuo(juo) - xuo(iuo)) 622 xc = xc + xij(1,ind) * recell(1,3) + & 623 xij(2,ind) * recell(2,3) + & 624 xij(3,ind) * recell(3,3) 625 if ( nint(xc) /= 0 ) cycle 626 end if 627 628 ! We also wish to remove the connection in 629 ! in the inner cell 630 ! I suspect we have a problem here 631 ! We may need a check for being in the unitcell 632 ! i.e.: 633 ! if ( l_RemUCellDistances ) then 634 ! if ( listh(ind) > no_u ) then ! As we check in the unit cell! 635 ! xo(1) = xa(1,liaorb(juo)) & 636 ! - xa(1,liaorb(iuo)) 637 ! xo(2) = xa(2,liaorb(juo)) & 638 ! - xa(2,liaorb(iuo)) 639 ! xo(3) = xa(3,liaorb(juo)) & 640 ! - xa(3,liaorb(iuo)) 641 ! else 642 ! xo = 0.0_dp 643 ! end if 644 ! end if 645 if ( l_RemUCellDistances ) then 646 xo(1) = xa(1,liaorb(juo)) & 647 - xa(1,liaorb(iuo)) 648 xo(2) = xa(2,liaorb(juo)) & 649 - xa(2,liaorb(iuo)) 650 xo(3) = xa(3,liaorb(juo)) & 651 - xa(3,liaorb(iuo)) 652 end if 653 654 kxij = & 655 k(1) * (xij(1,ind) - xo(1)) + & 656 k(2) * (xij(2,ind) - xo(2)) + & 657 k(3) * (xij(3,ind) - xo(3)) 658 cphase = exp(dcmplx(0d0,kxij)) 659 Hk(iuo,juo) = Hk(iuo,juo)+H(ind)*cphase 660 Sk(iuo,juo) = Sk(iuo,juo)+S(ind)*cphase 661 end do 662 end do 663 664 else setup_HS 665 666 do iuo = 1,no_tot 667 iu = iuo + l_RemNFirstOrbitals 668 iind = listhptr(iu) 669 do j = 1,numh(iu) 670 ind = iind + j 671 juo = listh(ind) - l_RemNFirstOrbitals 672 673 ! Cycle if we are not in the middle region 674 if ( juo < 1 .or. no_tot < juo ) cycle 675 676 ! If we have a removal of the z-direction do the following 677 if ( l_RemZConnection ) then 678 xc = - (xuo(juo) - xuo(iuo)) 679 xc = xc + xij(1,ind) * recell(1,3) + & 680 xij(2,ind) * recell(2,3) + & 681 xij(3,ind) * recell(3,3) 682 if ( nint(xc) /= 0 ) cycle 683 end if 684 685 Hk(iuo,juo) = Hk(iuo,juo)+H(ind) 686 Sk(iuo,juo) = Sk(iuo,juo)+S(ind) 687 end do 688 end do 689 690 end if setup_HS 691 692 if ( l_RemUCellDistances ) then 693 call memory('D','I',no_tot,'set_HS') 694 deallocate(liaorb) 695 end if 696 697 if ( l_RemZConnection ) then 698 call memory('D','D',no_tot,'set_HS') 699 deallocate(xuo) 700 end if 701 702 end subroutine set_HS_matrix_2d 703 704!***************** 705! Setting the Hamiltonian for a specific k-point as a specific transfer 706! This can be used to obtain the transfer matrix for a certain Hamilton and S. 707! It requires that the Node has the full sparse matrix available 708!***************** 709 subroutine set_HS_transfermatrix_1d(Gamma,ucell,na_u,no_u,no_s,maxnh, & 710 xij,numh,listhptr,listh,H,S, & 711 k,transfer_cell,HkT,SkT,xa,iaorb, & 712 DUMMY, & ! Ensures that the programmer makes EXPLICIT keywork passing 713 RemUCellDistances,RemNFirstOrbitals,RemNLastOrbitals) 714 use sys, only : die 715 use alloc, only : re_alloc 716 use geom_helper, only : ucorb 717 use cellSubs, only : reclat 718 719! *********************** 720! * INPUT variables * 721! *********************** 722 logical, intent(in) :: Gamma ! Is it a Gamma Calculation? 723 real(dp), intent(in) :: ucell(3,3) ! The unit cell of system 724 integer, intent(in) :: na_u ! Unit cell atoms 725 integer, intent(in) :: no_u ! Unit cell orbitals 726 integer, intent(in) :: no_s ! Supercell orbitals 727 integer, intent(in) :: maxnh ! Hamiltonian size 728 real(dp), intent(in) :: xij(3,maxnh) ! differences with unitcell, differences with unitcell 729 integer, intent(in) :: numh(no_u),listhptr(no_u) 730 integer, intent(in) :: listh(maxnh) 731 real(dp), intent(in) :: H(maxnh) ! Hamiltonian 732 real(dp), intent(in) :: S(maxnh) ! Overlap 733 real(dp), intent(in) :: k(3) ! k-point in [1/Bohr] 734 integer, intent(in) :: transfer_cell(3) ! The transfer cell directions 735 real(dp), intent(in) :: xa(3,na_u) ! Atomic coordinates (needed for RemZConnection & RemUCellDistances) 736 integer, intent(in) :: iaorb(no_u) ! The equivalent atomic index for a given orbital (needed for RemUCellDistances) 737! *********************** 738! * OUTPUT variables * 739! *********************** 740 complex(dp), pointer :: HkT(:), SkT(:) 741 742! *********************** 743! * OPTIONAL variables * 744! *********************** 745 logical, intent(in), optional :: DUMMY ! Do not supply this, it merely requires the coder 746! ! to use the keyworded arguments! 747 logical, intent(in), optional :: RemUCellDistances 748 integer, intent(in), optional :: RemNFirstOrbitals, RemNLastOrbitals 749 750! *********************** 751! * LOCAL variables * 752! *********************** 753 real(dp) :: recell(3,3) ! The reciprocal unit cell 754 real(dp) :: xo(3), xijo(3), xc 755 real(dp) :: kxij 756 complex(dp) :: cphase 757 integer :: no_tot 758 integer :: i,j,iu,iuo,juo,iind,ind 759 logical :: l_RemUCellDistances 760 integer :: l_RemNFirstOrbitals, l_RemNLastOrbitals 761 762 if ( present(DUMMY) ) & 763 call die("You must specify the keyworded arguments & 764 &for set_HS_transfermatrix") 765 766 ! Option collecting 767 l_RemUCellDistances = .false. 768 if ( present(RemUCellDistances) ) & 769 l_RemUCellDistances = RemUCellDistances 770 771 ! Make l_RemNFirstOrbitals contain the number of orbitals 772 ! to be removed from the Hamiltonian in the BEGINNING 773 l_RemNFirstOrbitals = 0 774 if ( present(RemNFirstOrbitals) ) & 775 l_RemNFirstOrbitals = RemNFirstOrbitals 776 ! Make l_RemNLastOrbitals contain the number of orbitals 777 ! to be removed from the Hamiltonian in the END 778 l_RemNLastOrbitals = 0 779 if ( present(RemNLastOrbitals) ) & 780 l_RemNLastOrbitals = RemNLastOrbitals 781 no_tot = no_u - (l_RemNLastOrbitals + l_RemNFirstOrbitals) 782 783 call re_alloc(HkT,1,no_tot*no_tot,name='HkT',routine='set_HS') 784 call re_alloc(SkT,1,no_tot*no_tot,name='SkT',routine='set_HS') 785 786 ! Prepare the cell to calculate the index of the atom 787 call reclat(ucell,recell,0) ! Without 2*Pi 788 789! 790! Setup H,S for this transfer k-point: 791! 792 do i = 1,no_tot*no_tot 793 HkT(i) = dcmplx(0.d0,0.d0) 794 SkT(i) = dcmplx(0.d0,0.d0) 795 end do 796 797 xo(:) = 0.0_dp 798 799 setup_HST: if (.not.Gamma ) then 800 801 do iuo = 1 , no_tot 802 iu = iuo + l_RemNFirstOrbitals 803 iind = listhptr(iu) 804 do j = 1,numh(iu) 805 ind = iind + j 806 juo = ucorb(listh(ind),no_u) - l_RemNFirstOrbitals 807 808 ! Cycle if we are not in the middle region 809 if ( juo < 1 .or. no_tot < juo ) cycle 810 811 ! Determine the transfer matrix cell 812 xijo(:) = xij(:,ind) - ( & 813 xa(:,iaorb(juo+l_RemNFirstOrbitals)) - xa(:,iaorb(iu)) & 814 ) 815 xc = sum(xijo(:) * recell(:,1)) 816 if ( nint(xc) /= transfer_cell(1) .and. transfer_cell(1) /= TRANSFER_ALL ) cycle 817 xc = sum(xijo(:) * recell(:,2)) 818 if ( nint(xc) /= transfer_cell(2) .and. transfer_cell(2) /= TRANSFER_ALL ) cycle 819 xc = sum(xijo(:) * recell(:,3)) 820 if ( nint(xc) /= transfer_cell(3) .and. transfer_cell(3) /= TRANSFER_ALL ) cycle 821 822 ! We also wish to remove the connection in 823 ! in the inner cell 824 ! I suspect we have a problem here 825 ! We may need a check for being in the unitcell 826 ! i.e.: 827 ! if ( l_RemUCellDistances ) then 828 ! if ( listh(ind) > no_u ) then ! As we check in the unit cell! 829 ! xo(1) = xa(1,iaorb(juo + l_RemNFirstOrbitals)) & 830 ! - xa(1,iaorb(iu)) 831 ! xo(2) = xa(2,iaorb(juo + l_RemNFirstOrbitals)) & 832 ! - xa(2,iaorb(iu)) 833 ! xo(3) = xa(3,iaorb(juo + l_RemNFirstOrbitals)) & 834 ! - xa(3,iaorb(iu)) 835 ! else 836 ! xo = 0.0_dp 837 ! end if 838 ! end if 839 if ( l_RemUCellDistances ) then 840 xo(1) = xa(1,iaorb(juo + l_RemNFirstOrbitals)) & 841 - xa(1,iaorb(iu)) 842 xo(2) = xa(2,iaorb(juo + l_RemNFirstOrbitals)) & 843 - xa(2,iaorb(iu)) 844 xo(3) = xa(3,iaorb(juo + l_RemNFirstOrbitals)) & 845 - xa(3,iaorb(iu)) 846 end if 847 848 kxij = & 849 k(1) * (xij(1,ind) - xo(1)) + & 850 k(2) * (xij(2,ind) - xo(2)) + & 851 k(3) * (xij(3,ind) - xo(3)) 852 cphase = exp(dcmplx(0d0,kxij)) 853 i = iuo+(juo-1)*no_tot 854 HkT(i) = HkT(i)+H(ind)*cphase 855 SkT(i) = SkT(i)+S(ind)*cphase 856 end do 857 end do 858 859 else setup_HST 860 ! It is not a Gamma calculation, thus we do not have any 861 ! neighbouring cells etc. 862 do iuo = 1 , no_tot 863 iu = iuo + l_RemNFirstOrbitals 864 iind = listhptr(iu) 865 do j = 1,numh(iu) 866 ind = iind + j 867 juo = listh(ind) - l_RemNFirstOrbitals 868 869 ! Cycle if we are not in the middle region 870 if ( juo < 1 .or. no_tot < juo ) cycle 871 872 ! Determine the transfer matrix cell 873 xijo(:) = xij(:,ind) - ( & 874 xa(:,iaorb(juo+l_RemNFirstOrbitals)) - xa(:,iaorb(iu)) & 875 ) 876 xc = sum(xijo(:) * recell(:,1)) 877 if ( nint(xc) /= transfer_cell(1) .and. transfer_cell(1) /= TRANSFER_ALL ) cycle 878 xc = sum(xijo(:) * recell(:,2)) 879 if ( nint(xc) /= transfer_cell(2) .and. transfer_cell(2) /= TRANSFER_ALL ) cycle 880 xc = sum(xijo(:) * recell(:,3)) 881 if ( nint(xc) /= transfer_cell(3) .and. transfer_cell(3) /= TRANSFER_ALL ) cycle 882 883 i = iuo+(juo-1)*no_tot 884 HkT(i) = HkT(i)+H(ind) 885 SkT(i) = SkT(i)+S(ind) 886 end do 887 end do 888 889 end if setup_HST 890 891 end subroutine set_HS_transfermatrix_1d 892 893 894 subroutine set_HS_transfermatrix_2d(Gamma,ucell,na_u,no_u,no_s,maxnh, & 895 xij,numh,listhptr,listh,H,S, & 896 k,transfer_cell,HkT,SkT,xa,iaorb,& 897 DUMMY, & ! Ensures that the programmer makes EXPLICIT keywork passing 898 RemUCellDistances,RemNFirstOrbitals,RemNLastOrbitals) 899 use sys, only : die 900 use alloc, only : re_alloc 901 use geom_helper, only : ucorb 902 use cellSubs, only : reclat 903 904! *********************** 905! * INPUT variables * 906! *********************** 907 logical, intent(in) :: Gamma ! Is it a Gamma Calculation? 908 real(dp), intent(in) :: ucell(3,3) ! The unit cell of system 909 integer, intent(in) :: na_u ! Unit cell atoms 910 integer, intent(in) :: no_u ! Unit cell orbitals 911 integer, intent(in) :: no_s ! Total orbitals 912 integer, intent(in) :: maxnh ! Hamiltonian size 913 real(dp), intent(in) :: xij(3,maxnh) ! differences with unitcell, differences with unitcell 914 integer, intent(in) :: numh(no_u),listhptr(no_u) 915 integer, intent(in) :: listh(maxnh) 916 real(dp), intent(in) :: H(maxnh) ! Hamiltonian 917 real(dp), intent(in) :: S(maxnh) ! Overlap 918 real(dp), intent(in) :: k(3) ! k-point in [1/Bohr] 919 integer, intent(in) :: transfer_cell(3) ! The transfer cell directions 920 real(dp), intent(in) :: xa(3,na_u) ! Atomic coordinates (needed for RemZConnection & RemUCellDistances) 921 integer, intent(in) :: iaorb(no_u) ! The equivalent atomic index for a given orbital (needed for RemUCellDistances) 922! *********************** 923! * OUTPUT variables * 924! *********************** 925 complex(dp), pointer :: HkT(:,:), SkT(:,:) 926! *********************** 927! * OPTIONAL variables * 928! *********************** 929 logical, intent(in), optional :: DUMMY ! Do not supply this, it merely requires the coder 930! ! to use the keyworded arguments! 931 logical, intent(in), optional :: RemUCellDistances 932 integer, intent(in), optional :: RemNFirstOrbitals, RemNLastOrbitals 933 934! *********************** 935! * LOCAL variables * 936! *********************** 937 real(dp) :: recell(3,3) 938 real(dp) :: xo(3), xijo(3), xc 939 integer :: no_tot 940 real(dp) :: kxij 941 complex(dp) :: cphase 942 integer :: j,iuo,iu,juo,iind,ind 943 logical :: l_RemUCellDistances 944 integer :: l_RemNFirstOrbitals, l_RemNLastOrbitals 945 946 if ( present(DUMMY) ) & 947 call die("You must specify the keyworded arguments & 948 &for set_HS_transfermatrix") 949 950 ! Option collecting 951 l_RemUCellDistances = .false. 952 if ( present(RemUCellDistances) ) & 953 l_RemUCellDistances = RemUCellDistances 954 955 ! Make l_RemNFirstOrbitals contain the number of orbitals 956 ! to be removed from the Hamiltonian in the BEGINNING 957 l_RemNFirstOrbitals = 0 958 if ( present(RemNFirstOrbitals) ) & 959 l_RemNFirstOrbitals = RemNFirstOrbitals 960 ! Make l_RemNLastOrbitals contain the number of orbitals 961 ! to be removed from the Hamiltonian in the END 962 l_RemNLastOrbitals = 0 963 if ( present(RemNLastOrbitals) ) & 964 l_RemNLastOrbitals = RemNLastOrbitals 965 no_tot = no_u - (l_RemNLastOrbitals + l_RemNFirstOrbitals) 966 967 968 call re_alloc(HkT,1,no_tot,1,no_tot,name='HkT',routine='set_HS') 969 call re_alloc(SkT,1,no_tot,1,no_tot,name='SkT',routine='set_HS') 970 971 972 ! Prepare the cell to calculate the index of the atom 973 call reclat(ucell,recell,0) ! Without 2*Pi 974 975 ! Setup H,S for this k-point: 976 do juo = 1,no_tot 977 do iuo = 1,no_tot 978 HkT(iuo,juo) = dcmplx(0.d0,0.d0) 979 SkT(iuo,juo) = dcmplx(0.d0,0.d0) 980 end do 981 end do 982 983 xo(:) = 0.0_dp 984 985 setup_HST: if (.not.Gamma ) then 986 987 do iuo = 1,no_tot 988 iu = iuo + l_RemNFirstOrbitals 989 iind = listhptr(iu) 990 do j = 1,numh(iu) 991 ind = iind + j 992 juo = ucorb(listh(ind),no_u) - l_RemNFirstOrbitals 993 994 ! Cycle if we are not in the middle region 995 if ( juo < 1 .or. no_tot < juo ) cycle 996 997 ! Determine the transfer matrix cell 998 xijo(:) = xij(:,ind) - ( & 999 xa(:,iaorb(juo+l_RemNFirstOrbitals)) - xa(:,iaorb(iu)) & 1000 ) 1001 xc = sum(xijo(:) * recell(:,1)) 1002 if ( nint(xc) /= transfer_cell(1) .and. transfer_cell(1) /= TRANSFER_ALL ) cycle 1003 xc = sum(xijo(:) * recell(:,2)) 1004 if ( nint(xc) /= transfer_cell(2) .and. transfer_cell(2) /= TRANSFER_ALL ) cycle 1005 xc = sum(xijo(:) * recell(:,3)) 1006 if ( nint(xc) /= transfer_cell(3) .and. transfer_cell(3) /= TRANSFER_ALL ) cycle 1007 1008 ! We also wish to remove the connection in 1009 ! in the inner cell 1010 ! I suspect we have a problem here 1011 ! We may need a check for being in the unitcell 1012 ! i.e.: 1013 ! if ( l_RemUCellDistances ) then 1014 ! if ( listh(ind) > no_u ) then ! As we check in the unit cell! 1015 ! xo(1) = xa(1,iaorb(juo + l_RemNFirstOrbitals)) & 1016 ! - xa(1,iaorb(iu)) 1017 ! xo(2) = xa(2,iaorb(juo + l_RemNFirstOrbitals)) & 1018 ! - xa(2,iaorb(iu)) 1019 ! xo(3) = xa(3,iaorb(juo + l_RemNFirstOrbitals)) & 1020 ! - xa(3,iaorb(iu)) 1021 ! else 1022 ! xo = 0.0_dp 1023 ! end if 1024 ! end if 1025 if ( l_RemUCellDistances ) then 1026 xo(1) = xa(1,iaorb(juo + l_RemNFirstOrbitals)) & 1027 - xa(1,iaorb(iu)) 1028 xo(2) = xa(2,iaorb(juo + l_RemNFirstOrbitals)) & 1029 - xa(2,iaorb(iu)) 1030 xo(3) = xa(3,iaorb(juo + l_RemNFirstOrbitals)) & 1031 - xa(3,iaorb(iu)) 1032 end if 1033 1034 kxij = & 1035 k(1) * (xij(1,ind) - xo(1)) + & 1036 k(2) * (xij(2,ind) - xo(2)) + & 1037 k(3) * (xij(3,ind) - xo(3)) 1038 cphase = exp(dcmplx(0d0,kxij)) 1039 HkT(iuo,juo) = HkT(iuo,juo)+H(ind)*cphase 1040 SkT(iuo,juo) = SkT(iuo,juo)+S(ind)*cphase 1041 end do 1042 end do 1043 1044 else setup_HST 1045 1046 do iuo = 1,no_tot 1047 iu = iuo + l_RemNFirstOrbitals 1048 iind = listhptr(iu) 1049 do j = 1,numh(iu) 1050 ind = iind + j 1051 juo = listh(ind) - l_RemNFirstOrbitals 1052 1053 ! Cycle if we are not in the middle region 1054 if ( juo < 1 .or. no_tot < juo ) cycle 1055 1056 ! Determine the transfer matrix cell 1057 xijo(:) = xij(:,ind) - ( & 1058 xa(:,iaorb(juo+l_RemNFirstOrbitals)) - xa(:,iaorb(iu)) & 1059 ) 1060 xc = sum(xijo(:) * recell(:,1)) 1061 if ( nint(xc) /= transfer_cell(1) .and. transfer_cell(1) /= TRANSFER_ALL ) cycle 1062 xc = sum(xijo(:) * recell(:,2)) 1063 if ( nint(xc) /= transfer_cell(2) .and. transfer_cell(2) /= TRANSFER_ALL ) cycle 1064 xc = sum(xijo(:) * recell(:,3)) 1065 if ( nint(xc) /= transfer_cell(3) .and. transfer_cell(3) /= TRANSFER_ALL ) cycle 1066 1067 HkT(iuo,juo) = HkT(iuo,juo)+H(ind) 1068 SkT(iuo,juo) = SkT(iuo,juo)+S(ind) 1069 end do 1070 end do 1071 1072 end if setup_HST 1073 1074 end subroutine set_HS_transfermatrix_2d 1075 1076 subroutine set_HS_available_transfers(ucell,na_u,xa,lasto,no_u,maxnh, & 1077 xij,numh,listhptr,listh,transfer_cell) 1078 use geom_helper, only : ucorb, iaorb 1079 use cellSubs, only : reclat 1080 1081! *********************** 1082! * INPUT variables * 1083! *********************** 1084 real(dp), intent(in) :: ucell(3,3) ! The unit cell of system 1085 integer, intent(in) :: na_u ! Unit cell atoms 1086 real(dp), intent(in) :: xa(3,na_u) ! Atomic coordinates (needed for RemZConnection & RemUCellDistances) 1087 integer, intent(in) :: lasto(0:na_u) ! last orbital number of equivalent atom 1088 integer, intent(in) :: no_u ! Unit cell orbitals 1089 integer, intent(in) :: maxnh ! Hamiltonian size 1090 real(dp), intent(in) :: xij(3,maxnh) ! differences with unitcell, differences with unitcell 1091 integer, intent(in) :: numh(no_u),listhptr(no_u) 1092 integer, intent(in) :: listh(maxnh) 1093! *********************** 1094! * OUTPUT variables * 1095! *********************** 1096 integer, intent(out) :: transfer_cell(2,3) 1097 1098! *********************** 1099! * LOCAL variables * 1100! *********************** 1101 real(dp) :: recell(3,3) 1102 real(dp) :: xijo(3), xc 1103 integer :: ia, ja 1104 integer :: i,j,iuo,juo,ind 1105 1106 ! Initialize the transfer cell to: 1107 transfer_cell(:,:) = 0 1108 1109 ! Prepare the cell to calculate the index of the atom 1110 call reclat(ucell,recell,0) ! Without 2*Pi 1111 1112 do iuo = 1 , no_u 1113 ia = iaorb(iuo,lasto) 1114 do j = 1 , numh(iuo) 1115 ind = listhptr(iuo) + j 1116 juo = ucorb(listh(ind),no_u) 1117 ja = iaorb(juo,lasto) 1118 xijo(:) = xij(:,ind) - ( xa(:,ja)-xa(:,ia) ) 1119 ! Loop over directions 1120 do i = 1 , 3 1121 ! recell is already without 2*Pi 1122 xc = sum(xijo(:) * recell(:,i)) 1123 transfer_cell(1,i) = min(transfer_cell(1,i),nint(xc)) 1124 transfer_cell(2,i) = max(transfer_cell(2,i),nint(xc)) 1125 end do 1126 end do 1127 end do 1128 1129 end subroutine set_HS_available_transfers 1130 1131 ! Routine for removing left right overlaps of certain regions. 1132 ! Is used to fully remove the connection between left and right 1133 ! states in the Hamiltonian 1134 subroutine matrix_rem_left_right_2d(no_tot,Hk,Sk,no_L,no_R) 1135 1136! ************************** 1137! * INPUT variables * 1138! ************************** 1139 integer, intent(in) :: no_tot, no_L, no_R 1140 1141! ************************** 1142! * OUTPUT variables * 1143! ************************** 1144 complex(dp), intent(inout) :: Hk(no_tot,no_tot), Sk(no_tot,no_tot) 1145 1146! ************************** 1147! * LOCAL variables * 1148! ************************** 1149 integer :: i,j 1150 1151 ! If nothing is to be removed return immidiately... 1152 if ( no_L == 0 .or. no_R == 0 ) return 1153 1154 do j = no_tot - no_R + 1 , no_tot 1155 do i = 1 , no_L 1156 Hk(i,j) = dcmplx(0.d0,0.d0) 1157 Sk(i,j) = dcmplx(0.d0,0.d0) 1158 Hk(j,i) = dcmplx(0.d0,0.d0) 1159 Sk(j,i) = dcmplx(0.d0,0.d0) 1160 end do 1161 end do 1162 1163 end subroutine matrix_rem_left_right_2d 1164 1165 subroutine matrix_rem_left_right_1d(no_tot,Hk,Sk,no_L,no_R) 1166 integer, intent(in) :: no_tot, no_L, no_R 1167 complex(dp), intent(inout) :: Hk(no_tot*no_tot), Sk(no_tot*no_tot) 1168 call matrix_rem_left_right_2d(no_tot,Hk,Sk,no_L,no_R) 1169 end subroutine matrix_rem_left_right_1d 1170 1171 1172 ! Routine for symmetrizing and shifting the matrix Ef 1173 subroutine matrix_symmetrize_2d(no_tot,Hk,Sk,Ef) 1174 1175! ************************** 1176! * INPUT variables * 1177! ************************** 1178 integer, intent(in) :: no_tot 1179 real(dp), intent(in) :: Ef 1180 1181! ************************** 1182! * OUTPUT variables * 1183! ************************** 1184 complex(dp), intent(inout) :: Hk(no_tot,no_tot), Sk(no_tot,no_tot) 1185 1186! ************************** 1187! * LOCAL variables * 1188! ************************** 1189 integer :: iuo,juo 1190 1191 do iuo = 1,no_tot 1192 do juo = 1,iuo-1 1193 1194 Sk(juo,iuo) = 0.5d0*( Sk(juo,iuo) + dconjg(Sk(iuo,juo)) ) 1195 Sk(iuo,juo) = dconjg(Sk(juo,iuo)) 1196 1197 Hk(juo,iuo) = 0.5d0*( Hk(juo,iuo) + dconjg(Hk(iuo,juo)) ) & 1198 - Ef*Sk(juo,iuo) 1199 Hk(iuo,juo) = dconjg(Hk(juo,iuo)) 1200 1201 end do 1202 1203 Sk(iuo,iuo)=Sk(iuo,iuo) - dcmplx(0d0,dimag(Sk(iuo,iuo)) ) 1204 1205 Hk(iuo,iuo)=Hk(iuo,iuo) - dcmplx(0d0,dimag(Hk(iuo,iuo)) ) & 1206 - Ef*Sk(iuo,iuo) 1207 end do 1208 1209 end subroutine matrix_symmetrize_2d 1210 1211 subroutine matrix_symmetrize_1d(no_tot,Hk,Sk,Ef) 1212 integer, intent(in) :: no_tot 1213 real(dp), intent(in) :: Ef 1214 complex(dp), intent(inout) :: Hk(no_tot*no_tot), Sk(no_tot*no_tot) 1215 call matrix_symmetrize_2d(no_tot,Hk(1),Sk(1),Ef) 1216 end subroutine matrix_symmetrize_1d 1217 1218end module m_hs_matrix 1219 1220