[Pw_forum] Re[Psi] and Im[Psi]

[O. Baris Malcioglu]

Dear Aritz,

First of all, you must notice that norm of the wavefunctions stored
will vary depending on what type of pseudopotential you use (Norm
conserving Ultrasoft etc.) The most straightforward in terms of
mathematical handling is the norm conserving one.

I prefer writing new subroutines rather than modifying old ones, since
I have witnessed doing that makes things get really messy after a
while.

Anyhow, here are some exempts from my personal notes, I hope you find
them useful. They are a bit unmaintained, thus as a disclaimer, please
report any error in them, so that I end up with better documentation.


1) How to calculate the real space grid points:

In this exempt, the aim is to find a set of grid points within a
certain distance from an atom.

The first step is to start an index rolling on available grid points

         DO ir = 1, nrxxs

nrxxs is explained as “the total dimension of the smooth grid” in
gsmooth (PW/pwcom.f90) . It is the total number of grid points in the
smooth grid.

            !
            ! ... three dimensional indexes
            !
            index = index0 + ir - 1
            k     = index / (nrx1s*nrx2s)
            index = index - (nrx1s*nrx2s)*k
            j     = index / nrx1s
            index = index - nrx1s*j
            i     = index
            !

Index0 is an offset for parallelism issues. It separates the grid into
processors and each processor starts processing from this offset.

nrx1s and nrx2s and nrx3s are again from gsmooth, and defined as
“maximum dimension of the smooth grid” in the corresponding direction,
they are used as the reduced forms of nrXs (X=1,2,3) dividing the grid
points among the nodes. In the serial case nrxxs=nrx1s*nrx2s*nrx3s but
this might or might not always be true for different parallel cases.
index (a temporary integer), nrx1s and nrx2s are all integers, so the
divisions above are integer divisions (otherwise the algorithm should
return 0 always, look at how the k and the index are calculated
consecutively).

As a convention, the segmentation of grid points should always start
from the 3rd dimension (z). The goal is to achieve an arrangement such
that the first 1..nrx1s elements correspond to ith element where j=0
k=0; the second 1..nrx1s element correspond to ith element where j=1
k=0 and etc.. This would make an index of higher rank increase only
when the cycles of all the lower ranks are complete. So the k index
(third dimension) is the total number of i and j cycles completed
(that would be nrx1s*nrx2s number of elements for each cycle). When k
times the number of elements per cycle is removed from the index
count, the excess only refers to i and j components. Using the same
trick to subtract the cycles of i (nrx1s elements each) the j index
value is obtained, and removing the number of elements per cycle for
each j increment from the index, the index for i is obtained.

            DO ipol = 1, 3
               posi(ipol) = DBLE( i )*inv_nr1s*at(ipol,1) + &
                            DBLE( j )*inv_nr2s*at(ipol,2) + &
                            DBLE( k )*inv_nr3s*at(ipol,3)
            END DO

this might be on a non-orthogonal basis and we should be extra careful
in the next step. Here ipol is an index running over the basis vectors
(these will be called x y and z directions for further reference)
at(ipol,1..3) are the simulation cell base vectors (in real space)
divided by alat (lattice parameter), in terms of Cartesian components
(second index). nr1s, nr2s and nr3s are the actual number of grid
points (where nrxXs are reduced in a node specific way to conserve
memory) thus at/nr is the vector increment for the particular
direction. Multiplying this with the increment counter in that
direction, a first iteration for the position vector is obtained.

            !
            posi(:) = posi(:) - tau_ia(:)

            !

This line offsets the position such that current atom is at the
center. tau_ia is just a rename for the tau from Modules/ions_base.f90
which contains the positions of the atomic cores.

            ! ... minimum image convenction
            !
            CALL cryst_to_cart( 1, posi, bg, -1 )
            !
            posi(:) = posi(:) - ANINT( posi(:) )
            !
            CALL cryst_to_cart( 1, posi, at, 1 )
            !

cryst_to_cart is a subroutine responsible for conversions of vectors
from cartesian coordinates to crystal coordinates (+1) and viceversa
(-1). This part ensures that any distance more than 0.5 crystal unit
far is treated according to the closest atomic center, not the one it
is further away. It is actually the minimum Brillouin cell for the
atom at hand.

            distsq = posi(1)**2 + posi(2)**2 + posi(3)**2
         END DO
      END DO

the distance from the center...

2) How to FFT a given band into real space (most simplistic version,
no gamma point tricks, no non-collinear, etc.)

            psic(nls(igk(:))) = evc(:,ibnd)
psic is the buffer where this kind of operations are performed.
Usually allocated in the initialization steps, so at most times can be
used directly. The igk(:) is the g,k correspondence array, and updated
for each k point under consideration. nls(igk(:)) returns the
corresponding index value in real space for the corresponding
reciprocal space point index value. evc are the wavefunctions. First
index is the grid point, second is the band.
             !
             CALL cft3s( psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2 )
internal pwscf wrapper for the fft that does the actual
transformation. After this, psic will contain the real space
wavefunctions with an index formalism for real space grids.


Best,
Baris Malcioglu
SISSA



On Wed, Mar 4, 2009 at 5:52 PM, Aritz Leonardo Liceranzu
<swblelia at ehu.es> wrote:
> Hi pwscf users,
>
> My question is regarding the output of the KS wave-functions.
>
> For a tddft problem that I am involved into, I need to numerically
> integrate the equilibrium KS wavefunctions. For this purpose I would
> like to obtain for each (k,i) Bloch state a chart containing:
>
> Rx,Ry,Rz, Re[psi], Im[psi]
>
> I surfed through old emails (like the one below) and I found that this
> feature is not implemented in PP codes (am I right?).
>
> Ab-init produces this kind of output but I am trying to avoid using it
> (I am loyal to espresso :) ). On the other hand, it shouldn't be very
> difficult to write a small program that does it, any hints so that I
> do not have to use ab-init????
>
> Thank you in advance
> Aritz
>
>
>
>
> Lorenzo Paulatto paulatto at sissa.it
> Sat Feb 23 12:48:15 CET 2008
>
> On Sab, 23 Febbraio 2008 5:42 am, nafise rezaei wrote:
>>> However, I want to extract and to plot parts of real and imaginary
>>> wave function.
>>> Do you know any program to generate wave function? If there is no
>>> program, I can modify what program to write it?
>
> The postprocessing utility pp.x can extract the square modulus of a
> wavefunction (you choose it by specifying kpoint and band index), plotting
> real and imaginary parts is not implemented, but you may be able to do it
> youself.
>
> On the other hand, if nobody have felt the need to implement it in the
> last not so few years there could be a reason. If you tell us what you
> want to do, instead of how you want to do it, you will probably get some
> more useful advice.
>
> Regards, LP
>
>
>
> _______________________________________________
> Pw_forum mailing list
> Pw_forum at pwscf.org
> http://www.democritos.it/mailman/listinfo/pw_forum
>