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 9 SUBROUTINE STM( NA, NO, NO_U, MAXNA, nspin,nspin_blocks,non_coll, 10 . ISA, IPHORB, INDXUO, LASTO, XA, CELL, UCELL, 11 . wf_unit, NK, gamma_wfsx, 12 . ZREF, ZMIN, ZMAX, NPX, NPY, NPZ, 13 . V0, EMAX, EMIN, 14 . ARMUNI, IUNITCD, RMAXO ) 15 16C ********************************************************************** 17C Simulate STM images in the Tersoff-Hamann approximation, by 18C extrapolating the wavefunctions into vacuum 19C 20C Coded by P. Ordejon and N. Lorente, November 2004 21C 22C Modified by N. Lorente, August 2005 23 ! Restructured by A. Garcia, March 2019 24C ********************************************************************** 25 26 use precision, only: dp, sp 27 USE ATMFUNCS 28 USE FDF 29 USE CHEMICAL 30 31 IMPLICIT NONE 32 33 real(dp), parameter :: Ang = 1.0_dp / 0.529177_dp 34 35 INTEGER, INTENT(IN) :: 36 . NA, NO, NO_U, NPX, NPY, NPZ, IUNITCD, 37 . nspin, nspin_blocks, MAXNA, NK, 38 . ISA(NA), IPHORB(NO), INDXUO(NO), LASTO(0:NA) 39 integer, intent(in) :: wf_unit 40 logical, intent(in) :: non_coll, gamma_wfsx 41 42 REAL(DP), INTENT(IN) :: 43 . ZMIN, ZMAX, ZREF, 44 . ARMUNI, RMAXO, V0, EMAX, EMIN 45 46 REAL(DP), INTENT(IN) :: CELL(3,3) 47 REAL(DP) :: UCELL(3,3), VOLCEL, XA(3,NA) 48 49 EXTERNAL :: VOLCEL 50 51C ****** INPUT ********************************************************* 52C INTEGER NA : Total number of atoms in Supercell 53C INTEGER NO : Total number of orbitals in Supercell 54C INTEGER NO_U : Total number of orbitals in Unit Cell 55C INTEGER MAXNA : Maximum number of neighbours of any atom 56C INTEGER ISA(NA) : Species index of each atom 57C INTEGER IPHORB(NO) : Orital index of each orbital in its atom 58C INTEGER INDXUO(NO) : Equivalent orbital in unit cell 59C INTEGER LASTO(0:NA) : Last orbital of each atom in array iphorb 60C REAL*8 XA(3,NA) : Atomic positions in cartesian coordinates 61C (in bohr), can be modified in the cube format 62C REAL*8 CELL(3,3) : Supercell vectors CELL(IXYZ,IVECT) 63C (in bohr) 64C REAL*8 UCELL(3,3) : Unit cell vectors CELL(IXYZ,IVECT) 65C (in bohr) 66C INTEGER NK : Number of k-points 67c REAL*8 ZREF : Position of reference plane for wf. estrapol. 68C REAL*8 ZMIN, ZMAX : Limits of the z-direction for the STM scan 69C INTEGER NPX,NPY,NPZ : Number of points along x and y and z 70C REAL*8 V0 : Value of the potential at the vacuum region in eV 71C REAL*8 EMAX : Maximum value for the energy window for STM in eV 72C REAL*8 EMIN : Minimum value for the energy window for STM in eV 73C REAL*8 ARMUNI : Conversion factor for the charge density 74C INTEGER IUNITCD : Unit of the charge density 75C REAL*8 RMAXO : Maximum range of basis orbitals 76C ********************************************************************** 77 78 INTEGER, DIMENSION(:), ALLOCATABLE :: JNA 79 REAL(DP), DIMENSION(:), ALLOCATABLE :: R2IJ 80 REAL(DP), DIMENSION(:,:), ALLOCATABLE :: XIJ 81 82 INTEGER 83 . IA, ISEL, NNA, I, J, IN, IAT1, IO, IUO, IAVEC1, 84 . IS1, IPHI1, NX, NY, NZ, IWF, IK, ISPIN, grid_u, str_u, 85 . IX, IY, IZ, NSX, NSY, NAU, iv, is 86 87 REAL(DP) 88 . DOT, RMAX, XPO(3), RMAX2, XVEC1(3), 89 . PHIMU, GRPHIMU(3), 90 . PHASE, SI, CO, ENER, PMIKR, SIMIKR, COMIKR, USAVE, VC, VU 91 92 real(dp) :: total_weight, k(3) 93 94 REAL(DP), ALLOCATABLE :: RHO(:,:,:,:) 95 REAL(SP), ALLOCATABLE :: wf_single(:,:) 96 COMPLEX(DP), ALLOCATABLE :: wf(:,:) 97 REAL(DP), ALLOCATABLE :: wk(:) 98 99 COMPLEX(DP) EXPPHI, EXMIKR, d11, d12, d21, d22 100 101 ! The last dimension of these is the number of spinor components 102 ! 1 for collinear, and 2 for NC/SOC 103 COMPLEX(DP), ALLOCATABLE :: CW(:,:,:), CWE(:,:,:,:), CWAVE(:) 104 105 LOGICAL FIRST 106 integer :: idummy, number_of_wfns, spinor_comps 107 integer :: nspin_rho ! Number of components needed for rho array ("spin%Grid") 108 109 CHARACTER :: SNAME*40, FNAME*256, stm_label*60 110 111 EXTERNAL :: NEIGHB, IO_ASSIGN, IO_CLOSE 112 113C ********************************************************************** 114C INTEGER IA : Atom whose neighbours are needed. 115C A routine initialization must be done by 116C a first call with IA = 0 117C If IA0=0, point X0 is used as origin instead 118C INTEGER ISEL : Single-counting switch (0=No, 1=Yes). If ISEL=1, 119C only neighbours with JA.LE.IA are included in JNA 120C INTEGER NNA : Number of non-zero orbitals at a point in 121C real space 122C INTEGER JNA(MAXNA) : Atom index of neighbours. The neighbours 123C atoms might be in the supercell 124C REAL*8 XIJ(3,MAXNA) : Vectors from point in real space to orbitals 125C REAL*8 R2IJ(MAXNA) : Squared distance to atomic orbitals 126C REAL*8 XPO(3) : Coordinates of the point of the plane respect 127C we are going to calculate the neighbours orbitals 128C INTEGER IZA(NA) : Atomic number of each atom 129C ********************************************************************** 130 131 132 ! The first dimension of wf_single is the number of real numbers per orbital 133 ! to be read from the WFSX file: 134 ! 1 for real wfs, 2 for complex, and four for the two spinor components 135 ! wf is a complex array which holds either a wfn or a two-component spinor. 136 137 if (non_coll) then 138 allocate(wf_single(4,1:no_u)) 139 allocate(wf(1:no_u,2)) 140 spinor_comps = 2 141 else 142 spinor_comps = 1 143 if (gamma_wfsx) then 144 allocate(wf_single(1,1:no_u)) 145 allocate(wf(1:no_u,1)) 146 else 147 allocate(wf_single(2,1:no_u)) 148 allocate(wf(1:no_u,1)) 149 endif 150 endif 151 152C Initialize neighbour subroutine -------------------------------------- 153 IA = 0 154 ISEL = 0 155 RMAX = RMAXO 156 NNA = MAXNA 157 158 IF (ALLOCATED(JNA)) THEN 159 DEALLOCATE(JNA) 160 ENDIF 161 IF (ALLOCATED(R2IJ)) THEN 162 DEALLOCATE(R2IJ) 163 ENDIF 164 IF (ALLOCATED(XIJ)) THEN 165 DEALLOCATE(XIJ) 166 ENDIF 167 168 ALLOCATE(JNA(MAXNA)) 169 ALLOCATE(R2IJ(MAXNA)) 170 ALLOCATE(XIJ(3,MAXNA)) 171 172 nspin_rho = min(4,nspin) 173 174 allocate(CWAVE(spinor_comps)) 175 ALLOCATE(CW(0:NPX-1,0:NPY-1,spinor_comps)) 176 ALLOCATE(CWE(0:NPX-1,0:NPY-1,0:NPZ-1,spinor_comps)) 177 ALLOCATE(RHO(0:NPX-1,0:NPY-1,0:NPZ-1,nspin_rho)) 178 179 FIRST = .TRUE. 180 DO I = 1,3 181 XPO(I) = 0.D0 182 ENDDO 183 CALL NEIGHB( CELL, RMAX, NA, XA, XPO, IA, ISEL, 184 . NNA, JNA, XIJ, R2IJ, FIRST ) 185 FIRST = .FALSE. 186 RMAX2 = RMAXO**2 187 188 IF (.not. monoclinic(ucell)) then 189 WRITE(6,*) 'error: the code only accepts monoclinic cells' 190 WRITE(6,*) ' with Z as the vertical axis' 191 STOP 192 ENDIF 193 194! Initialize density 195 196 RHO = 0 197 198! Loop over k-points and wavefunctions to include in the STM image 199 200 allocate(wk(nk)) 201 DO IK = 1, NK 202 do ispin = 1, nspin_blocks 203 read(wf_unit) idummy, k(1:3), wk(ik) 204 if (idummy /= ik) then 205 write(6,*) "ik index mismatch in WFS file" 206 WRITE(6,*) "ik in file, ik: ", idummy, ik 207 endif 208 read(wf_unit) idummy 209 if (idummy /= ispin) then 210 write(6,*) "ispin index mismatch in WFS file" 211 WRITE(6,*) "ispin in file, ispin: ", idummy, ispin 212 endif 213 read(wf_unit) number_of_wfns 214 215 WRITE(6,*) 'stm: Processing kpoint ',IK 216 WRITE(6,*) 'stm: nwf: ', number_of_wfns 217 WRITE(6,*) ' --------------------------------' 218 DO IWF = 1, number_of_wfns 219 read(wf_unit) idummy 220 if (idummy /= iwf) then 221 ! The file holds a subset of wfs, with the original indexes... 222 WRITE(6,*) 'Original wf index: ', idummy 223 endif 224 read(wf_unit) ener 225 226 ! Check that we have a bound state (E below vacuum level), 227 ! in the chosen window 228 229 IF (ENER .LT. EMIN .OR. ENER .GT. EMAX) then 230 read(wf_unit) ! skip wfn info 231 CYCLE 232 ENDIF 233 234 IF (ENER .GT. V0) THEN 235 WRITE(6,*) 'ERROR: ENERGY EIGENVALUE ',IWF, 236 . ' FOR K-POINT ', IK, 'FOR SPIN ',ISPIN 237 WRITE(6,*) ' IS ABOVE VACUUM LEVEL' 238 STOP 239 ENDIF 240 241 WRITE(6,"(a,i5,i2)") 'stm: wf (spin) in window: ',iwf,ispin 242 243 read(wf_unit) (wf_single(:,io), io=1,no_u) 244 ! Use a double precision complex form in what follows 245 if ( non_coll) then 246 wf(:,1) = cmplx(wf_single(1,:), wf_single(2,:), kind=dp) 247 wf(:,2) = cmplx(wf_single(3,:), wf_single(4,:), kind=dp) 248 else 249 if (gamma_wfsx) then 250 wf(:,1) = cmplx(wf_single(1,:), 0.0_sp, kind=dp) 251 else 252 wf(:,1) = cmplx(wf_single(1,:),wf_single(2,:),kind=dp) 253 endif 254 endif 255 256 ! Loop over all points in real space 257 ! The last point (zmax) is now included by making stepz=(Zmax-Zmin)/(NPZ-1) 258 ! This forces the definition of a slightly larger c vector below. 259 ! In this way, the last plane recorded in the file will correspond to Z=Zmax 260 DO NZ = 0, NPZ-1 261 262 if (npz == 1) then 263 XPO(3) = ZMIN 264 else 265 XPO(3) = ZMIN + NZ*(ZMAX-ZMIN)/(NPZ-1) 266 endif 267 268 if ( XPO(3) < Zref ) then 269 ! Initialize density to unextrapolated density 270 271 WRITE(6,"(a,f10.4)") 'stm: Using plain LDOS for z =', 272 $ xpo(3) 273 DO NY = 0,NPY-1 274 DO NX = 0,NPX-1 275 276 ! Note that the (periodic) X and Y directions 277 ! are treated as usual, with smaller step 278 XPO(1) = NX*UCELL(1,1)/NPX + 279 $ NY*UCELL(1,2)/NPY 280 XPO(2) = NX*UCELL(2,1)/NPX + 281 $ NY*UCELL(2,2)/NPY 282 283 call get_cwave(wf(:,1:spinor_comps)) 284 285 ! Now for the various cases 286 if (nspin <= 2) then 287 RHO(NX,NY,NZ,ispin) = RHO (NX,NY,NZ,ispin) 288 & + REAL(CWAVE(1)*CONJG(CWAVE(1)), dp) 289 $ * ARMUNI * WK(IK) 290 else ! non-collinear 291 ! CHECK THIS 292 d11 = cwave(1) * conjg(cwave(1)) 293 d12 = cwave(1) * conjg(cwave(2)) 294 d21 = cwave(2) * conjg(cwave(1)) 295 d22 = cwave(2) * conjg(cwave(2)) 296 297 ! Hermitify? 298 D12 = 0.5_dp * (D12 + conjg(D21)) 299 300 ! Recall: dm(:,3) = real(d12); dm(:,4) = -aimag(d12) 301 rho(nx,ny,nz,1) = rho(nx,ny,nz,1) 302 $ + real(d11,dp) * armuni * wk(ik) 303 rho(nx,ny,nz,2) = rho(nx,ny,nz,2) 304 $ + real(d22,dp) * armuni * wk(ik) 305 rho(nx,ny,nz,3) = rho(nx,ny,nz,3) 306 $ + real(d12,dp) * armuni * wk(ik) 307 rho(nx,ny,nz,4) = rho(nx,ny,nz,4) 308 $ - aimag(d12) * armuni * wk(ik) 309 endif 310 ENDDO 311 ENDDO 312 313 else 314 315 ! Extrapolate from reference plane 316 ! Compute value of the wfn at this reference plane 317 WRITE(6,"(a,i4)") 'stm: Extrapolating from nz:', nz 318 319 DO NY = 0,NPY-1 320 DO NX = 0,NPX-1 321 322 XPO(1) = NX*UCELL(1,1)/NPX + 323 $ NY*UCELL(1,2)/NPY 324 XPO(2) = NX*UCELL(2,1)/NPX + 325 $ NY*UCELL(2,2)/NPY 326 XPO(3) = ZREF 327 328 call get_cwave(wf(:,1:spinor_comps)) 329 CW(NX,NY,1:spinor_comps) = 330 $ CWAVE(1:spinor_comps) 331 332 ENDDO 333 ENDDO 334 335 ! This is mildly wasteful in terms of initialization 336 ! of FFTW, but it will do for now 337 ! We assume that both spinor components are propagated 338 ! in the same way 339 do is = 1, spinor_comps 340 CALL EXTRAPOLATE(NPX,NPY,NPZ,ZREF,ZMIN,ZMAX, 341 $ UCELL,V0, 342 . CW(0,0,is),ENER,K,CWE(0,0,0,is)) 343 enddo 344 345 ! Now for the various cases 346 ! Be careful not to overwrite the z<zref parts... 347 if (nspin <= 2) then 348 RHO(:,:,NZ:,ispin) = RHO (:,:,NZ:,ispin) 349 & + REAL(CWE(:,:,NZ:,1)* 350 $ CONJG(CWE(:,:,NZ:,1)), dp) 351 $ * ARMUNI * WK(IK) 352 else ! non-collinear 353 354 do iz = nz, npz-1 355 do ny=0,npy-1 356 do nx=0,npx-1 357 d11 = cwe(nx,ny,iz,1) * conjg(cwe(nx,ny,iz,1)) 358 d12 = cwe(nx,ny,iz,1) * conjg(cwe(nx,ny,iz,2)) 359 d21 = cwe(nx,ny,iz,2) * conjg(cwe(nx,ny,iz,1)) 360 d22 = cwe(nx,ny,iz,2) * conjg(cwe(nx,ny,iz,2)) 361 362 ! Hermitify? 363 D12 = 0.5_dp * (D12 + conjg(D21)) 364 365 ! Recall: dm(:,3) = real(d12); dm(:,4) = -aimag(d12) 366 rho(nx,ny,iz,1) = rho(nx,ny,iz,1) 367 $ + real(d11,dp) * armuni * wk(ik) 368 rho(nx,ny,iz,2) = rho(nx,ny,iz,2) 369 $ + real(d22,dp) * armuni * wk(ik) 370 rho(nx,ny,iz,3) = rho(nx,ny,iz,3) 371 $ + real(d12,dp) * armuni * wk(ik) 372 rho(nx,ny,iz,4) = rho(nx,ny,iz,4) 373 $ - aimag(d12) * armuni * wk(ik) 374 enddo 375 enddo 376 enddo 377 endif ! collinear or not 378 379 ! And we are done with the z planes 380 EXIT ! loop over NZ 381 382 endif ! z below or above Zref 383 384 ENDDO ! NZ 385 ENDDO ! ispin = 1, nspin_blocks 386 387 ENDDO ! wfn number 388 ENDDO ! k-point 389 390 ! This should not be necessary if a proper BZ-sampled set of wfs is used 391 total_weight = sum(wk(1:nk)) 392 rho = rho / total_weight 393 ! Normalize if not spin-polarized 394 if ((nspin_blocks == 1) .and. (.not. non_coll)) then 395 rho = 2.0_dp * rho 396 endif 397 398! Write charge density in Siesta format 399 400 call io_assign(grid_u) 401 SNAME = FDF_STRING('SystemLabel','siesta') 402 stm_label = FDF_STRING('stm-label','') 403 if (stm_label == '') then 404 FNAME = trim(SNAME) // '.STM.LDOS' 405 else 406 FNAME = trim(SNAME) // '.' // trim(stm_label) // '.STM.LDOS' 407 endif 408 409 WRITE(6,*) 410 WRITE(6,*) 'stm: writing SIESTA format file ', FNAME 411 WRITE(6,*) 412 413 open(grid_u,file=FNAME,form='unformatted', 414 . status='unknown') 415 ! write the correct range of the z axis: temporarily override ucell 416 USAVE = UCELL(3,3) 417 ! Make the cell slightly taller, so that the last (NPZ-1) plane corresponds 418 ! to Z=Zmax 419 if (npz == 1 ) then 420 UCELL(3,3) = 1.0_dp ! We have a single plane. Use a 1.0 bohr-thick height 421 else 422 UCELL(3,3) = NPZ * abs(ZMAX-ZMIN) / (NPZ-1) 423 endif 424 425 WRITE(grid_u) UCELL, 426 $ [0.0_dp, 0.0_dp, ZMIN], ! Extra info for origin 427 $ [.true.,.true.,.false.] ! Periodic ? 428 429 WRITE(grid_u) NPX, NPY, NPZ, nspin_rho 430 431 do ispin = 1, nspin_rho 432 DO IZ=0,NPZ-1 433 DO IY=0,NPY-1 434 WRITE(grid_u) (REAL(RHO(IX,IY,IZ,ispin),sp),IX=0,NPX-1) 435 ENDDO 436 ENDDO 437 enddo 438 call io_close(grid_u) 439 440 ! Write a dummy STRUCT file for use with the 'grid to cube' converter 441 call io_assign(str_u) 442 open(str_u,file=trim(SNAME)//'.CELL_STRUCT', 443 $ form='formatted', status='unknown') 444 write(str_u,'(3x,3f18.9)') ((UCELL(ix,iv)/Ang,ix=1,3),iv=1,3) 445 write(str_u,*) 0 ! number of atoms 446 close(str_u) 447 448 ! restore ucell 449 UCELL(3,3)=USAVE 450 451 DEALLOCATE(RHO) 452 DEALLOCATE(CWE) 453 DEALLOCATE(CW) 454 deallocate(cwave) 455 456 CONTAINS 457 458 subroutine get_cwave(psi) 459 complex(dp), intent(in) :: psi(:,:) ! Can deal with spinors 460 461 ! Inherits all data by host association 462 463! The periodic part of the Bloch functions is defined by 464! \begin{equation} 465! u_{n \vec{k}} (\vec{r}) = 466! \sum_{\vec{R} \mu} c_{n \mu}(\vec{k}) 467! e^{i \vec{k} \cdot ( \vec{r}_{\mu} + \vec{R} - \vec{r} )} 468! \phi_{\mu} (\vec{r} - \vec{r}_{\mu} - \vec{R} ) , 469!\end{equation} 470! 471!\noindent where $\phi_{\mu} (\vec{r} - \vec{r}_{\mu} - \vec{R} )$ 472! is an atomic orbital of the basis set centered on atom $\mu$ in 473! the unit cell $\vec{R}$, and $c_{n \mu}(\vec{k})$ are the coefficients 474! of the wave function 475 476 CWAVE = (0.0D0, 0.0D0) 477 478! CWAVE is meant to be the periodic part of the wavefunction, 479! for both the computation of the charge density directly (exp(ikr) phase 480! is irrelevant) and for propagation of the wave function (?) 481 482 ! First step (see below) 483 ! Phase to cancel the phase of the wave function: -i.k.r 484 485 PMIKR = -(K(1)*XPO(1) + K(2)*XPO(2) + K(3)*XPO(3)) 486 SIMIKR=DSIN(PMIKR) 487 COMIKR=DCOS(PMIKR) 488 EXMIKR=DCMPLX(COMIKR,SIMIKR) 489 490C Localize non-zero orbitals at each point in real space --------------- 491 492 IA = 0 493 ISEL = 0 494 NNA = MAXNA 495 ! Get neighbors of point xpo 496 CALL NEIGHB( CELL, RMAX, NA, XA, XPO, IA, ISEL, 497 . NNA, JNA, XIJ, R2IJ, FIRST ) 498 499C Loop over Non-zero orbitals ------------------------------------------ 500 ! NOTE: If the z-extent of the box is large enough, we might be getting 501 ! contributions from orbitals in the next periodic image of the slab. 502 ! We should get rid of them 503 DO IAT1 = 1, NNA 504 IF( R2IJ(IAT1) .GT. RMAX2 ) CYCLE 505 506 IAVEC1 = JNA(IAT1) 507 IS1 = ISA(IAVEC1) 508 XVEC1(:) = -XIJ(:,IAT1) ! position of XPO with respect to 509 ! the atom 510 511 ! XPO + XIJ(IAT1) is just the absolute position of atom IAT1 512 ! We could cancel the phase above and keep only k*xij 513 514 PHASE = K(1)*(XPO(1)+XIJ(1,IAT1))+ 515 . K(2)*(XPO(2)+XIJ(2,IAT1))+ 516 . K(3)*(XPO(3)+XIJ(3,IAT1)) 517 518 SI=DSIN(PHASE) 519 CO=DCOS(PHASE) 520 EXPPHI=DCMPLX(CO,SI) 521 522 DO IO = LASTO(IAVEC1-1) + 1, LASTO(IAVEC1) 523 IPHI1 = IPHORB(IO) 524 IUO = INDXUO(IO) 525 CALL PHIATM( IS1, IPHI1, XVEC1, PHIMU, GRPHIMU ) 526 527 ! Note implicit loop over spinor components 528 CWAVE(:) = CWAVE(:) + PHIMU * psi(iuo,:) * EXPPHI * EXMIKR 529 530 ENDDO 531 ENDDO 532 533 end subroutine get_cwave 534 535 function monoclinic(cell) 536 real(dp), intent(in) :: cell(3,3) 537 logical monoclinic 538 539 real(dp), parameter :: tol = 1.0e-8_dp 540 541 ! This is too naive. It should be checking that a_3 is orthogonal 542 ! to both a_1 and a_2, but it is fine for this use case. 543 544 monoclinic = (abs(CELL(3,1)) < tol 545 $ .and. abs(CELL(3,2)) < tol 546 $ .and. abs(CELL(1,3)) < tol 547 $ .and. abs(CELL(2,3)) < tol ) 548 549 end function monoclinic 550 551 END 552