[Pw_forum] STS spectra

[Guido Fratesi]

Dear Daniele,

I'm sorry to reply only now to your message which dates two months ago, 
and I guess you could have probably found a convenient solution to get the 
LDOS by now.

Sometime ago I was also interested in computing the LDOS at fixed 
coordinate, r, and at variance with the energy, E.  I added a few bits to 
the code (3.1.1 at the time) which might be useful to you as well, but 
this never come to a stage sufficiently refined to be distributed and 
eventually got obsolete.

I attach here a few lines of the algorithm which I was adding into 
projwfc.f90, as you might be curious about that. Only the simplest, 
norm-conserving, collinear case was implemented. The basic idea is the 
following (possibly not the most efficient one), and closely follows the 
one for computing the projected DOS:

* given regions where to integrate the LDOS
* define "theta" arrays (dimension nrxx) which evaluate 1 inside and 0
  outside these regions --> thetathisproc(1:nrxx,ibox)
* for each band, each kpoint --> proj(ibox,ibnd,ik)
  - get the |wfc^2| (array of dimension nrxx)
  - multiply times the "theta" array and sum over space
  (this step is analogous to performing the
  projection over the atomic orbitals)
* sum over kpoints and bands --> ldosbox(ie,ibox,current_spin)

I would be interested to know if you find this way of proceeding 
convenient to your case (and in case I can send a sort-of-working example) 
or if you have spotted a better way to proceed, in the while.

Regards,
Guido

---------------------- compunting proj

  k_loop: DO ik = 1, nks
     !
     IF ( lsda ) current_spin = isk(ik)
     CALL gk_sort (xk (1, ik), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
     CALL davcio (evc, nwordwfc, iunwfc, ik, - 1)
     !
     bnd_loop: DO ibnd = 1, nbnd
        !
        IF (noncolin) THEN
           !
           CALL errore( 'projwfc', 'boxes and noncolin not implemented', 1 )
           !
        ELSE
           !
           caux(1:nrxx) = (0._DP,0._DP)
           DO ig = 1, npw
              caux (nl (igk (ig) ) ) = evc (ig, ibnd)
           ENDDO
           IF (gamma_only) THEN
              DO ig = 1, npw
                 caux (nlm(igk (ig) ) ) = CONJG(evc (ig, ibnd))
              ENDDO
           END IF
           CALL cft3 (caux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
           !
           raux(:) = DBLE( caux(:) )**2 + AIMAG( caux(:) )**2
           !
        END IF
        !
        ! The contribution of this wavefunction to the LDOS
        ! integrated in the volume is the projection of the
        ! squared wfc on a function =1 in the volume itself:
        !
        DO ibox = iboxmin, n_proj_boxes
           proj(ibox,ibnd,ik) = DDOT(nrxx,raux(:),1,thetathisproc(:,ibox),1) &
                &               / (nrx1*nrx2*nrx3)
        END DO
        !
     END DO bnd_loop
     !
     CALL reduce ((n_proj_boxes-iboxmin+1)*nbnd , proj(iboxmin,1,ik))
     !
  END DO k_loop

---------------------- summing over kp and bands

  DO ik = 1,nkstot
     IF ( nspin == 2 ) current_spin = isk ( ik )
     DO ibnd = 1, nbnd
        etev = et(ibnd,ik)
        ie_mid = NINT( (etev-Emin)/DeltaE )
        DO ie = MAX(ie_mid-ie_delta, 0), MIN(ie_mid+ie_delta, ne)
           delta = w0gauss((Emin+DeltaE*ie-etev)/degauss,ngauss) &
                 / degauss / rytoev
           !
           DO ibox = iboxmin, n_proj_boxes
              ldosbox(ie,ibox,current_spin)                = &
                   ldosbox(ie,ibox,current_spin)           + &
                   wk(ik) * delta * proj (ibox, ibnd, ik)
           END DO
           !
           ! dostot(:,ns) = total DOS (states/eV) for spin "ns"
           !
           dostot(ie,current_spin) = dostot(ie,current_spin) + &
                wk(ik) * delta
        END DO
     END DO
  END DO

----------------------



On Tue, 8 Jul 2008, Daniele Passerone wrote:

 | 
Dear Pwscf users and developers,

The postprocessing facilities of espresso allow to extract stm images, 
that correspond in Tersoff-Hamann approximation to the localized density 
of states integrated between Efermi and a given (positive or negative) 
bias. The more convenient way to analyze this quantity is in my opinion by 
generating a 3d grid and visualizing it with xcrysden, which incidentally 
can also show planes parallel to the cell sides.

The ILDOS option of pp.x (plot_num=10) somewhat extends this, allowing to 
compute the LDOS integrated between emin and emax. STS (scanning tunneling 
spectroscopy) experiments sample the local density of states at a surface: 
one obtains a profile of dI/dV at a certain point on the surface, and this 
gives valuable information about the local conductivity and presence of 
peaks in the DOS etc.

The question is: what is the more convenient and fast way to get the LDOS 
profile from Emin to Emax, on a certain plane (let us say, 2 Angstrom 
above the surface)? I thought one could make several jobs of pp.x : from 
Emin to Emin+DeltaE, from Emin + Delta E to Emin + 2*Delta E.... and so 
on... but maybe there are much better ways to do that.

Thank you in advance for any advice,
Daniele


-- 
Guido Fratesi

Dipartimento di Scienza dei Materiali
Universita` degli Studi di Milano-Bicocca