[Pw_forum] Correlated wave-function in real space

[Andreas Linscheid]

Dear Matteo,

Thank you for your quick respond!
My question is quite unusual, and I hope that I can clarify what I want.

We expand the superconducting order parameter
 chi = <psi^(r) psi^(r')>
in a Kohn-Sham orbital basis and solve the Kohn-Sham-Bogoliobov-deGennes 
equations for the expansion coefficients. What I doing now, is to 
multiply the basis with this expansion coefficients to recover the 
original quantity, chi. Therefore I am fine with the Kohn Sham states in 
their real space representation.
We want to have this quantity in a center of mass coordinate system R= 
r+r', S=r-r'. Since this is a 6 dim. quantity we fix S to a certain 
vector and then have something like the charge to plot.

I am writing a new module (routine) which computes this quantity chi for 
fixed S. This uses the module chdens.f90 to plot the desired quantity, 
like every other PW routine (I know).

Because of the not-working part of my program does essentially take the 
1D index of the wavefunction and knows what 3D vector this (should) 
correspond to,
ir => (i,j,k) => R
 then adds the 3D vector,
R => R +S,
 to obtain a new 1D index corresponding to R+S
R+S => ir'
(This is the crucial point), I was doing some experiments, as described 
below, that failed. Ok, the question is:

Is it correct that
1. The 3D real space vector R in lattice coordinates, that corresponds 
    to ir, where
    ir = k + (j-1)*nrx3s +(i-1)*nrx3s*nrx2s
    is
    R(1) = DBLE(i)/DBLE(nrx1s)
    R(2) = DBLE(j)/DBLE(nrx2s)
    R(3) = DBLE(k)/DBLE(nrx3s)
    ? (I am using nrx1s instead of nr1s on purpose!)

2. Is it correct, that the 3D vector to the same index ir in karthesian   
coordinates is
    R(:) = DBLE(i)/DBLE(nrx1s)*at(:,1) + DBLE(j)/DBLE(nrx2s)*at(:,2) + 
DBLE(k)/DBLE(nrx3s)*at(:,3)
? (I think in FAQ, this is described with nr3s,)

3. Can I take this vector (for example in lattice coords) and add
R(:) = R(:) + S(:)
map back to first unit cell
R(:) = R(:)-nint(R(:))
and map it back onto the real space FFT grid point (i_n,j_n,k_n):
i_n =nint(nrx1s*R(1)) + 1
j_n =nint(nrx2s*R(2)) + 1
k_n =nint(nrx3s*R(3)) + 1
What then has the the 1D index:
ir_n = k_n + (j_n-1)*nrx3s +(i_n-1)*nrx3s*nrx2s

Because if all this is correct I have the error somewhere else.
I hope I have clarified my problem?!
Thank you very much,

Andreas


Matteo Cococcioni wrote:
> Dear Andreas
>
> I'm not sure I understood correctly what you want to do. I have a couple 
> of questions/comments:
>
> 1) are you trying to reconstruct localized wavefunctions from the 
> Kohn-Sham states? This is what Wannier functions are all about.
>
> 2) If you want a translated wfc from the fourier transform of the 
> original one, you just have to multiply
> each Fourier transform at a given k by e^ik*S where S is the the amount 
> of translation you want to apply.
> However when you bring it back in real space I guess you will obtain 
> only the part of the translated wavefunction
> that is inside the first unit cell.
>
> Hope this helps.
>
> Matteo
>
>
> Andreas Linscheid wrote:
>   
>> 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)
>> _______________________________________________
>> Pw_forum mailing list
>> Pw_forum at pwscf.org
>> http://www.democritos.it/mailman/listinfo/pw_forum
>>   
>>     
>
>
>