1! 2! Copyright (C) 2003 A. Smogunov 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC). 9! 10! 11subroutine compbs(lleft, nocros, norb, nchan, kval, kfun, & 12 kfund, kint, kcoef, ikk, ien) 13! 14! Using the basis functions obtained by scatter_forw it computes 15! the complex band structure (CBS) of the lead. 16! Some variables needed for wave-function matching in transmission 17! calculation are constructed and saved. 18! 19 USE constants, ONLY : tpi 20 USE noncollin_module, ONLY : noncolin, npol 21 USE spin_orb, ONLY : lspinorb 22 USE lsda_mod, ONLY : nspin 23 USE cond 24 USE cell_base, ONLY : alat, at, omega 25 USE ions_base, ONLY : nat, ityp 26 27 implicit none 28 integer :: & 29 nocros, & ! number of orbitals crossing the boundary 30 noins, & ! number of interior orbitals 31 norb, & ! total number of orbitals 32 lleft ! 1/0 if it is left/right tip 33 integer :: ik, ikk, i, j, ig, n, iorb, iorb1, & 34 iorb2, aorb, borb, nchan, & 35 ij, is, js, ichan, ien 36 REAL(DP), PARAMETER :: eps=1.d-8 37 REAL(DP) :: raux, ddot 38 REAL(DP), ALLOCATABLE :: zpseu(:,:,:), zps(:,:) 39 COMPLEX(DP), PARAMETER :: cim=(0.d0,1.d0) 40 COMPLEX(DP) :: x1, & 41 kval(2*(n2d+npol*nocros)), kfun(n2d,2*(n2d+npol*nocros)), & 42 kint(nocros*npol,2*(n2d+npol*nocros)), & 43 kcoef(nocros*npol,2*(n2d+npol*nocros)), & 44 kfund(n2d,2*(n2d+npol*nocros)) 45 COMPLEX(DP), ALLOCATABLE :: amat(:,:), bmat(:,:), vec(:,:), & 46 zpseu_nc(:,:,:,:), zps_nc(:,:), & 47 aux(:,:), veceig(:,:), korb(:,:) 48 COMPLEX(DP), PARAMETER :: one=(1.d0,0.d0), zero=(0.d0,0.d0) 49 50 51 call start_clock('compbs') 52 53 noins = norb-2*nocros 54 IF (norb>0) THEN 55 if(lleft.eq.1) then 56 if (noncolin) then 57 allocate( zpseu_nc(2,norb,norb,nspin) ) 58 zpseu_nc = zpseul_nc 59 else 60 allocate( zpseu(2,norb,norb) ) 61 zpseu = zpseul 62 endif 63 else 64 if (noncolin) then 65 allocate( zpseu_nc(2,norb,norb,nspin) ) 66 zpseu_nc = zpseur_nc 67 else 68 allocate( zpseu(2,norb,norb) ) 69 zpseu = zpseur 70 endif 71 endif 72 if (noncolin) then 73 allocate( zps_nc( norb*npol, norb*npol ) ) 74 else 75 allocate( zps( norb, norb ) ) 76 endif 77 END IF 78 79 allocate( amat( (2*n2d+npol*norb), (2*n2d+npol*norb) ) ) 80 allocate( bmat( (2*n2d+npol*norb), (2*n2d+npol*norb) ) ) 81 allocate( vec( (2*n2d+npol*norb), 2*(n2d+npol*nocros) ) ) 82 allocate( aux( n2d, 2*n2d+npol*norb)) 83 IF (lorb) allocate( korb(npol*(nocros+noins),2*(n2d+npol*nocros)) ) 84 85 amat=(0.d0,0.d0) 86 bmat=(0.d0,0.d0) 87 88! 89! zps=zpseu-e*qq for US-PP and zps=zpseu for norm-conserv. PP 90! 91 do iorb=1, norb 92 do iorb1=1, norb 93 if (noncolin) then 94 ij=0 95 do is=1,npol 96 do js=1,npol 97 ij=ij+1 98 zps_nc(npol*(iorb-1)+is, npol*(iorb1-1)+js)= & 99 zpseu_nc(1,iorb,iorb1,ij) 100 if (lspinorb) then 101 zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= & 102 zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) & 103 -eryd*zpseu_nc(2,iorb,iorb1,ij) 104 else 105 if (is.eq.js) & 106 zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= & 107 zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) & 108 -eryd*zpseu_nc(2,iorb,iorb1,ij) 109 endif 110 enddo 111 enddo 112 else 113 zps(iorb,iorb1)=zpseu(1,iorb,iorb1)-eryd*zpseu(2,iorb,iorb1) 114 endif 115 enddo 116 enddo 117 118! 119! Forming the matrices A and B for generalized eigenvalue problem 120 121! 122! 1 123! 124 do n=1, 2*n2d 125 do ig=1, n2d 126 amat(ig, n)=fun1(ig, n) 127 amat(ig+n2d,n)=fund1(ig, n) 128 bmat(ig,n)=fun0(ig, n) 129 bmat(ig+n2d,n)=fund0(ig, n) 130 enddo 131 enddo 132 133! 134! 2 135! 136 do iorb=1, norb*npol 137 do ig=1, n2d 138 amat(ig, 2*n2d+iorb)=funl1(ig, iorb) 139 amat(n2d+ig, 2*n2d+iorb)=fundl1(ig, iorb) 140 bmat(ig, 2*n2d+iorb)=funl0(ig, iorb) 141 bmat(n2d+ig, 2*n2d+iorb)=fundl0(ig, iorb) 142 enddo 143 enddo 144! 145! 3 146! 147 do iorb=1, norb*npol 148 aorb=iorb 149 borb=iorb 150 if (iorb.le.npol*nocros) aorb=iorb+npol*(noins+nocros) 151 if (iorb.gt.npol*nocros) borb=iorb-npol*(noins+nocros) 152 do n=1, 2*n2d 153 do iorb1=1, norb*npol 154 if (noncolin) then 155 amat(2*n2d+iorb,n)=amat(2*n2d+iorb,n)+ & 156 zps_nc(aorb, iorb1)*intw1(iorb1,n) 157 if (borb.gt.0) bmat(2*n2d+iorb,n)= & 158 bmat(2*n2d+iorb,n)-zps_nc(borb,iorb1)*intw1(iorb1,n) 159 else 160 amat(2*n2d+iorb,n)=amat(2*n2d+iorb,n)+ & 161 zps(aorb,iorb1)*intw1(iorb1,n) 162 if (borb.gt.0) bmat(2*n2d+iorb,n)= & 163 bmat(2*n2d+iorb,n)-zps(borb,iorb1)*intw1(iorb1,n) 164 endif 165 enddo 166 enddo 167 enddo 168! 169! 4 170! 171 do iorb=1, nocros*npol 172 do iorb1=1, norb*npol 173 do iorb2=1, norb*npol 174 if (noncolin) then 175 bmat(2*n2d+iorb,2*n2d+iorb1)=bmat(2*n2d+iorb,2*n2d+iorb1) & 176 -zps_nc(iorb,iorb2)*intw2(iorb2, iorb1) 177 else 178 bmat(2*n2d+iorb,2*n2d+iorb1)=bmat(2*n2d+iorb,2*n2d+iorb1) & 179 -zps(iorb,iorb2)*intw2(iorb2, iorb1) 180 endif 181 enddo 182 bmat(2*n2d+iorb+npol*(noins+nocros),2*n2d+iorb1)= & 183 bmat(2*n2d+iorb,2*n2d+iorb1) 184 enddo 185 bmat(2*n2d+iorb,2*n2d+iorb)= & 186 bmat(2*n2d+iorb,2*n2d+iorb)+(1.d0,0.d0) 187 enddo 188! 189! 5 190! 191 do iorb=1, norb*npol 192 aorb=iorb 193 if (iorb.le.npol*nocros) aorb=iorb+npol*(noins+nocros) 194 do iorb1=1, norb*npol 195 do iorb2=1, norb*npol 196 if (noncolin) then 197 amat(2*n2d+iorb,2*n2d+iorb1)=amat(2*n2d+iorb,2*n2d+iorb1)+ & 198 zps_nc(aorb,iorb2)*intw2(iorb2, iorb1) 199 else 200 amat(2*n2d+iorb,2*n2d+iorb1)=amat(2*n2d+iorb,2*n2d+iorb1)+ & 201 zps(aorb,iorb2)*intw2(iorb2, iorb1) 202 endif 203 enddo 204 enddo 205 if (aorb.eq.iorb) amat(2*n2d+iorb,2*n2d+iorb)= & 206 amat(2*n2d+iorb,2*n2d+iorb)-(1.d0,0.d0) 207 enddo 208 209! 210! To reduce matrices and solve GEP A X = c B X; X = {a_n, a_\alpha} 211! 212 213 call compbs_2(npol*nocros, npol*norb, n2d, 2*(n2d+npol*nocros), & 214 amat, bmat, vec, kval) 215 216! 217! To normalize (over XY plane) all the states 218! 219 call zgemm('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb, & 220 one, amat, 2*n2d+npol*norb, vec, 2*n2d+npol*norb, & 221 zero, kfun, n2d) 222 do ig=1,n2d 223 do ik=1, 2*n2d+npol*norb 224 aux(ig,ik)=amat(n2d+ig,ik) 225 enddo 226 enddo 227 call zgemm('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb, & 228 one, aux, n2d, vec, 2*n2d+npol*norb, zero, kfund, n2d) 229 230 do ik=1, 2*(n2d+npol*nocros) 231 raux=ddot(2*n2d,kfun(1,ik),1,kfun(1,ik),1)*sarea 232 raux=1.d0/sqrt(raux) 233 call dscal(2*(2*n2d+npol*norb),raux,vec(1,ik),1) 234 call dscal(2*n2d,raux,kfun(1,ik),1) 235 call dscal(2*n2d,raux,kfund(1,ik),1) 236 enddo 237 238! 239! To find k-vector and the current of Bloch states 240! 241 242 call kbloch (2*(n2d+npol*nocros), kval) 243 244 call jbloch(2*(n2d+npol*nocros), n2d, norbf, norb, nocros, & 245 kfun, kfund, vec, kval, intw1, intw2, nchan, npol) 246! 247! To save band structure result 248! 249 kfun=(0.d0,0.d0) 250 kfund=(0.d0,0.d0) 251 kint=(0.d0,0.d0) 252! 253! To account for the case of the right lead 254! 255 if (lleft.eq.0) then 256 do i=1, 2*n2d 257 do j=1, 2*n2d+npol*norb 258 amat(i,j)=bmat(i,j) 259 enddo 260 enddo 261 do i=2*n2d+1, 2*n2d+npol*nocros 262 do j=1, 2*n2d+npol*norb 263 amat(i,j)=-bmat(i+npol*(nocros+noins),j) 264 enddo 265 enddo 266 endif 267! 268! psi_k and psi'_k on the scattering region boundary 269! 270 call zgemm('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb,& 271 one, amat, 2*n2d+npol*norb, vec, 2*n2d+npol*norb, & 272 zero, kfun, n2d) 273 do ig=1,n2d 274 do ik=1, 2*n2d+npol*norb 275 aux(ig,ik)=amat(n2d+ig,ik) 276 enddo 277 enddo 278 call zgemm('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb,& 279 one, aux, n2d, vec, 2*n2d+npol*norb, zero, kfund, n2d) 280 281! 282! kint(iorb, ik)=\sum_{iorb1} D_{iorb,iorb1} 283! \int_{cell} W_iorb1^* psi_ik ) 284! for the orbitals crossing the boundary 285! 286 do ik=1, 2*(n2d+npol*nocros) 287 do iorb=1, nocros*npol 288 do j=1, 2*n2d+npol*norb 289 kint(iorb,ik)=kint(iorb,ik)+amat(2*n2d+iorb,j)*vec(j,ik) 290 enddo 291 enddo 292 enddo 293! 294! a_iorb = kcoef(iorb,ik) = \sum_{iorb1} D_{iorb,iorb1} 295! \int_{all space} W_iorb1^* psi_ik ) 296! for the orbitals crossing the boundary 297! 298 do ik=1, 2*(n2d+npol*nocros) 299 do iorb=1, nocros*npol 300 if (lleft.eq.0) then 301 kcoef(iorb,ik)=vec(2*n2d+iorb,ik) 302 else 303 kcoef(iorb,ik)=vec(2*n2d+npol*(nocros+noins)+iorb,ik) 304 endif 305 enddo 306 enddo 307 308! 309! to set up B.S. for the right lead in the case of identical tips 310! 311 if(lleft.eq.1.and.ikind.eq.1) then 312 nchanr=nchan 313 call dcopy(2*(n2d+npol*nocros), kval, 1, kvalr, 1) 314 kfunr=(0.d0,0.d0) 315 kfundr=(0.d0,0.d0) 316 kintr=(0.d0,0.d0) 317 318 do i=1, 2*n2d 319 do j=1, 2*n2d+npol*norb 320 amat(i,j)=bmat(i,j) 321 enddo 322 enddo 323 do i=2*n2d+1, 2*n2d+npol*nocros 324 do j=1, 2*n2d+npol*norb 325 amat(i,j)=-bmat(i+npol*(nocros+noins),j) 326 enddo 327 enddo 328 329 do ik=1, 2*(n2d+npol*nocros) 330 do ig=1, n2d 331 do j=1, 2*n2d+npol*norb 332 kfunr(ig,ik)= kfunr(ig,ik)+amat(ig,j)*vec(j,ik) 333 kfundr(ig,ik)= kfundr(ig,ik)+amat(n2d+ig,j)*vec(j,ik) 334 enddo 335 enddo 336 enddo 337 do ik=1, 2*(n2d+npol*nocros) 338 do iorb=1, nocros*npol 339 do j=1, 2*n2d+npol*norb 340 kintr(iorb,ik)=kintr(iorb,ik)+amat(2*n2d+iorb,j)*vec(j,ik) 341 enddo 342 enddo 343 enddo 344 do ik=1, 2*(n2d+npol*nocros) 345 do iorb=1, nocros*npol 346 kcoefr(iorb,ik)=vec(2*n2d+iorb,ik) 347 enddo 348 enddo 349 350 endif 351 352!-- 353! integrals of Bloch states with boundary orbitals for left/right leads 354! 355 if (lorb) then 356 korb = 0.d0 357 do ik = 1, 2*(n2d+npol*nocros) 358 359 do iorb = 1, npol*nocros 360 iorb1 = iorb + npol*(nocros+noins) 361 do ig = 1, 2*n2d 362 korb(iorb,ik) = korb(iorb,ik)+ & 363 intw1(iorb1,ig)*vec(ig,ik) 364 enddo 365 do ig = 1, norb*npol 366 korb(iorb,ik) = korb(iorb,ik)+ & 367 intw2(iorb1,ig)*vec(2*n2d+ig,ik) 368 enddo 369 enddo 370 371 do iorb = 1, npol*nocros 372 x1 = 0.d0 373 do ig = 1, 2*n2d 374 x1 = x1 + intw1(iorb,ig)*vec(ig,ik) 375 enddo 376 do ig = 1, norb*npol 377 x1 = x1 + intw2(iorb,ig)*vec(2*n2d+ig,ik) 378 enddo 379 korb(iorb,ik) = korb(iorb,ik) + x1* & 380 exp(kval(ik)*(0.d0,1.d0)*tpi) 381 enddo 382 enddo 383 384 if (ikind.ne.2.or.lleft.ne.0) korbl(:,:) = korb(:,:) 385 if (ikind.ne.2.or.lleft.ne.1) then 386 do ik = 1, 2*(n2d+npol*nocros) 387 x1 = exp(-kval(ik)*(0.d0,1.d0)*tpi) 388 do iorb = 1, npol*nocros 389 korbr(iorb,ik) = korb(iorb,ik) * x1 390 enddo 391 enddo 392 endif 393 394 endif 395!-- 396 397!-- 398! Computes and writes the propagating Bloch states 399! 400 if (lorb.and.ikind.eq.0.and.nchan /= 0) then 401 allocate( veceig(nchan, nchan) ) 402 deallocate( aux ) 403 allocate( aux(4*n2d+npol*(norb+2*nocros),nchan) ) 404 405!-- right moving states 406 veceig = 0.d0 407 aux = 0.d0 408 do ichan = 1, nchan 409 do ig = 1, 2*n2d + npol*norb 410 aux(ig, ichan) = vec(ig,ichan) 411 enddo 412 enddo 413 CALL scat_states_plot(ikk,ien,norb,nocros,nchan,aux,veceig,.true.) 414 415!-- left moving states 416 veceig = 0.d0 417 aux = 0.d0 418 do ichan = 1, nchan 419 do ig = 1, 2*n2d + npol*norb 420 aux(ig, ichan) = vec(ig,n2d+npol*nocros+ichan) 421 enddo 422 enddo 423 CALL scat_states_plot(ikk,ien,norb,nocros,nchan,aux,veceig,.false.) 424 425 deallocate( veceig ) 426 endif 427!-- 428 429 deallocate(amat) 430 deallocate(bmat) 431 deallocate(vec) 432 deallocate(aux) 433 IF (norb>0) THEN 434 if (noncolin) then 435 deallocate(zpseu_nc) 436 deallocate(zps_nc) 437 else 438 deallocate(zpseu) 439 deallocate(zps) 440 endif 441 ENDIF 442 if (lorb) deallocate(korb) 443 call stop_clock('compbs') 444 445 return 446end subroutine compbs 447