[Pw_forum] Correlated wave-function in real space

[Andreas Linscheid]

Dear PWSCF users,

I am currently writing an extension to the post processing part, 
combining information by an other program, to compute the 
superconducting order parameter in real space. In order to do this I 
need to calculate a correlated wave-function in real space, like

CONJG(psi_nk(R - Constvect)) * psi_nk(R + Constvect),

where R and Constvect are 3D real space vectors and Constvect is held 
fixed. I didn't succeed up to now, so I ask for help.
My final goal is to find the corresponding FFT index, ir <=> R so that I 
can add (or subtract) S: R -> R+(-)S <=> ir_plus(minus)
I figured out that the error in my program is where I touch the index 
ir, so it seems I misunderstood something of the correspondence between 
1D consequtively ordered grid index ir and 3D real space vector R. My 
underlying question is therefore:

How can I evaluate a given wavefunction psi(ir,ibnd) in real space at a 
given 3D vector R?

What I tried up to now:
1. To get a better understanding of the wavefunction in real space, I 
tried to reproduce the sum of wavefunction (absolute value squared),

Result(ir) = sum(irreducible grid k,bands) coeff_nk*|wavefunc_at_k(ir)|^2

using the normal ouput (chdens.f90), by avoiding the 
fourier-interpolation in chdens.f90 and plotting it directly. (This means:)
I compute the wavefunction like in local_dos.f90.

call gk_sort (xk(1,ik), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
!load the PW coefficients of the hit irreducible k point
call davcio (evc, nwordwfc, iunwfc, kfull_to_irred(ik), - 1)

do ig = 1, npw
    do ibnd = 1, nbnd
        wavefunc_at_k ( nls(igk(ig)),ibnd) = evc (ig, ibnd)
    enddo
enddo

call cft3s (wavefunc_at_k(:,ibnd), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)

I get the right result, i.e. I can reproduce the charge density when 
setting all coefficients above the fermi level to zero, only when I 
don't touch the index ir.

2. According to the FAQ I assumed that after the FFT, wavefunc_at_k ( 
ir,ibnd) is the wavefunction of band iband at the cartesian coordinate

R(:) = DBLE(i-1)/DBLE(nrx1s)*at(:,1) + DBLE(j-1)/DBLE(nrx2s)*at(:,2) + 
DBLE(k-1)/DBLE(nrx3s)*at(:,3)
when
ir = k + (j-1)*nrx3s +(i-1)*nrx3s*nrx2s

Therefore I though that I should be able to reproduce the  Result above 
(maybe in a different numerical resolution) when simply writing directly 
a file.

k = nrx3s/2 !look at the boron plane
do i = 1, nrx1s
  do j = 1, nrx2s
    Rvect(:) = DBLE(i-1)/DBLE(nrx1s)*at(:,1) + 
DBLE(j-1)/DBLE(nrx2s)*at(:,2) + DBLE(k-1)/DBLE(nrx3s)*at(:,3)
      ir = 1 + (j-1)*nrx3s + (i-1)*nrx3s*nrx2s + (k-1)
      write(13,'(2f10.5,f15.8)') Rvect(1),Rvect(2),Result(ir)
      write(13,*)
  enddo
enddo

Well I get something different which does not even have the symmetry 
group of the lattice. I also tried this in lattice coords:

R(:) = DBLE(i-1)/DBLE(nrx1s)*at(:,1) + DBLE(j-1)/DBLE(nrx2s)*at(:,2) + 
DBLE(k-1)/DBLE(nrx3s)*at(:,3)
call cryst_to_cart(1,R,bg,-1)

Neither did this work (of cause the same problem). If anyone can 
confirm, that my thoughts are (in principle) right, I at least know that 
the bug must be somewhere else.

I very much appreciate your help!!

Andreas Linscheid
(Fu-Berlin, AG-Gross)