[Pw_forum] possible bug in scatter_forw.f90

[Alexander]

Dear Manoj

On Tuesday 21 October 2008 00:01, Manoj Srivastava wrote:
> Dear Alexander,
>  Thank you very much for your answer. I have one more question, it seems
> my list of questions is never gonna end :(
It is always a pleasure when there is someone going through
the code written by you and trying to understand what is there 
after all :-).

>  this is about init_gper.f90 subroutine. where you calculate
>
> |G+k|^2<Ecut2d.
>
> the expression -
> norm2=(((il+xyk(1,ik))*bg(1,1)+(jl+xyk(2,ik))*bg(1,2))**2+ &
>             ((il+xyk(1,ik))*bg(2,1)+(jl+xyk(2,ik))*bg(2,2))**2)* &
>
> i believe is  |G+k|^2. How do you get above expression in the code for
> this? 
norm2 is indeed |G_perp+k_perp|^2, where G_perp=il*b_1+jl*b_2 and 2d vector 
k_perp=xyk(1,ik)*b_1+xyk(2,ik)*b_2 and you indeed select out those 
G_perp+k_perp which satisfy (G_perp+k_perp)^2<Ecut2d. 
The number of latters is ngper. Notice that gper array contains 
the cartesian coordinates of these G_perp+k_perp vectors.
> I tried to work it out but could not get this. You even have g.k 
> term!  
> I guess i would have more understanding once i understand how did 
> you define Gper in the fourier transform,
>  V(rper,z)=sum(gper)[exp(i*gper*rper) V(gper,z)
> Is gper just gx and gy, and then you rotate your axis to get
> bg(1,1),bg(1,2), bg(2,1) and bg(2,2)? Would you mind explaining?
The potential V(z,G_perp) is calculated for the full set of 2d G_perp related 
to the real space discretization of r_perp=(x,y). 
V(z,G_perp) <--FFT--> V(z,r_perp).
G_perp=il*b_1+jl*b_2, il=1..nrx, jl=1..nry, where nrx and nry are in fact
nr1s and nr2s of pw.x and are defined by the cut-off ecutwfc.
This potential is then used in the subroutine local.f90.

Hope this helps,
Alexander



>
> Regards,
> Manoj
>
> On Mon, 20 Oct 2008, Alexander wrote:
> > Dear Manoj
> >
> > On Thursday 16 October 2008 23:15, Manoj Srivastava wrote:
> > > Dear Alexandar,
> > >  Thank you very much. Your answers make sense and now i have better
> > > understanding of the subroutine scatter_forw.f90.  I have one question
> > > about poten.f90 subroutine, where you calculate the V(p,gper) in each
> > > slices of the slab(equation 15). I have a question about Gz, which is
> > > 2*pi*q/L, with [-N/2]+1<=q>=[N/2]. Now in the code at the line at the
> > > line 60, you have Gz multiplied with bg(3,3). I kind of understand this
> > > is because 3rd direction is perpendicular to xy plane,
> >
> > Exactly. As it is implemented now pwcond deals with systems
> > having monoclinic unit cell, i.e. such that a3 is perpendicular to both
> > a1 and a2 (a1 and a2 are not necessary orthogonal one to another) and the
> > direction of transport is a3. This seems to us the most general setup for
> > transport and I do not think we should consider at the moment the
> > situation where the a3 is not perpendicular to the xy plane.
> >
> > Best regards,
> > Alexander
> >
> > > but not totally clear
> > > about it. Would you mind explaining? my other question is when
> > > reciprocal space vectors b1, b2 and b3 are not orhogonal to each other,
> > > which is more general case, then should not we take bg(1,3) and bg(2,3)
> > > in there too?
> > >  Thank you once again for listening to my questions and i hope i am
> > > not bothering you too much.
> > >
> > > Regards,
> > > Manoj Srivastava
> > > Physics Graduate Student
> > > Department of Physics
> > > Gainesville, FL, USA.
> > >
> > > have  On Wed, 15 Oct 2008, Alexander wrote:
> > > > Dear Manoj
> > > >
> > > > On Tuesday 14 October 2008 17:49, Manoj Srivastava wrote:
> > > > > Dear Alexader,
> > > > >
> > > > >  Thank you very much for your reply. I was out of town and could
> > > > > not see the message till yesterday. I have some questions on it, it
> > > > > would be great if you could answer them.
> > > > >
> > > > > On Fri, 10 Oct 2008, Alexander wrote:
> > > > > > Dear Manoj
> > > > > >
> > > > > > On Friday 10 October 2008 05:20, Manoj Srivastava wrote:
> > > > > > > Any follow ups on the message below? Alexander?
> > > > > > >
> > > > > > > Regards,
> > > > > > > Manoj
> > > > > > >
> > > > > > > On Wed, 8 Oct 2008, Manoj Srivastava wrote:
> > > > > > > > Dear All,
> > > > > > > >  Is any body familiar with the scatter_forw.f90 subroutine of
> > > > > > > > PWCOND? I think there is a bug in this subroutine at the
> > > > > > > > place where you calculate intw2, which is z integration of
> > > > > > > > nonlocal wavefunction with beta function. I have looked up
> > > > > > > > the necessary formulae in the Choi & Ihm's paper (PRB 59,
> > > > > > > > 2267, Jan 1999). In the paper, nonlocal wavefunction has 3
> > > > > > > > terms, one of which contains beta function, say W. The other
> > > > > > > > two parts dont have any W in them, rather they are just plane
> > > > > > > > wave solution. (Please have a look at equation 24 and 26 of
> > > > > > > > the paper ). So, when we are doing z integration of nonlocal
> > > > > > > > wavefunction with beta function W, we should have three
> > > > > > > > terms, one of which should contain two W, but rest should
> > > > > > > > just have one W. On the other hand in the code all three
> > > > > > > > terms contain two beta functions!! (please have a look at
> > > > > > > > line 228 of scatter_forw.f90). I am wondering if my
> > > > > > > > understanding is right?
> > > > > >
> > > > > > Your understanding is wrong. On the line 228 you add to the
> > > > > > integral ONLY contribution with two beta functions (from the 1st
> > > > > > term in the nonlocal function, see Eq. 24) but on the line 393
> > > > > > the contribution from the 2nd term is added too. There is no 3rd
> > > > > > term since every step you rotate your local solutions \psi_n in
> > > > > > such a way that b_{\lambda \alpha lm} vanish (see the paragraph
> > > > > > after Eq. 37).
> > > > >
> > > > >   I believe that the equations 40, 41 and 42 implemented in the
> > > > > code are in the momentum space. In the code, intw2 has 3 terms. We
> > > > > can tell this by looking at the subroutine integrals.f90, where
> > > > > int2d is defined. Now, by looking at the choi& Ihm 's paper, the
> > > > > very first term of psi{alpha,l,m} has
> > > > > f{lambda,alpha,l,m}(equation26).And, if we substitute this in
> > > > > equation 40, we should just get one term. I am just confuesd on the
> > > > > fact that how come this one term in paper gets transformed into 3
> > > > > terms of the code!
> > > >
> > > > intw2(alpha,beta) is the integral of nonlocal solution (24) with w
> > > > functions, intw2(alpha,beta) = \int_{0,d}
> > > > w_{alpha}(r)*psi_{beta}(r)dr, where alpha,beta={alpha,l,m} in the
> > > > paper. In each slab it will have TWO contributions from the first and
> > > > second terms of (24), the one from the third term vanishes as I
> > > > explained before. These two contributions are calculated on the lines
> > > >  228  and 393 of the code. The integrals intw1 (alpha,n)=\int_{0,d}
> > > > w_{alpha}(r)*psi_{n}(r)dr are the integrals of local solutions (17)
> > > > with the same projector functions w.
> > > > The integrals intw2 and intw1 have nothing to do with boundary
> > > > conditions (40). They enter in (11) and are needed to construct the
> > > > final solutions according to Eq.(10). Substituting the form of this
> > > > solution (10) in the boundary conditions (40) you will be led to the
> > > > generalized eigenvalue problem, AX=exp(ikd) BX, with known matrices A
> > > > and B and unknown coefficients X={a_n,C_{alpha,l,m}} defining the
> > > > resulting solution (10).
> > > >
> > > > > Would you mind explaining? Also, the integrals that are defined in
> > > > > the code are only in the slabs (nz1), so how are you taking the
> > > > > integration in the region (0,d), which we need according to
> > > > > equations 40, 41 and 42 ?
> > > >
> > > > When you go over slabs you add new contributions to the integrals
> > > > intw1 and intw2 (see the line 228 for the first contribution to intw2
> > > > and the part "add to the integrals" for another contribution to intw2
> > > > and contributions to intw1) so at the end you obtain integrals over
> > > > all the region.
> > > >
> > > > Regards,
> > > > Alexander
> > > >
> > > > > >  > > I have one more question in the later part of the same
> > > > > >  > > subroutine. What
> > > > > > > >
> > > > > > > > does the lapack subroutine
> > > > > > > > ZGESV(2*n2d,2*n2d+norb*npol,amat,2*n2d,ipiv,xmat,2*n2d,info)
> > > > > > > > do?
> > > > > > > > I tried to look it up and i got the idea that it is trying to
> > > > > > > > solve amat*x=xmat, with amat and xmat known and x unknown,
> > > > > > > > and at the end of calculation it stores x in xmat.
> > > > > > > > so, basically it does--  x=[(amat)^(-1)]*xmat. Am i right?
> > > > > > > > So, it changes the structure of xmat from what is defined in
> > > > > > > > 'constructs matrices' part. Is it correct?
> > > > > >
> > > > > > Here, you are right.
> > > > > >
> > > > > > > > Also, afterwards in this code where  it 'rotates integrals'
> > > > > > > > is not very clear to me.
> > > > > > > > Could somebody please tell me in little detail, what is going
> > > > > > > > on here?
> > > > > >
> > > > > > As I explained before rotation means that every step (at each
> > > > > > slab) you make a linear combination of local solutions \psi_n
> > > > > > with the transformation matrix  h_{nn'} (see the last paragraph
> > > > > > on the p. 2270). so that the new solutions have the property:
> > > > > > b_{\lambda n} = \delta_{\lambda n} and the new nonlocal solutions
> > > > > > have b_{\lambda \alpha lm}=0. Doing this transformation of local
> > > > > > and nonlocal solutions you should also perform the corresponding
> > > > > > transformations of the integrals.
> > > > > >
> > > > > > > > Also, is this subroutine written just on the basis of Choi
> > > > > > > > and Ihm paper, or are there more reference to it? If yes,
> > > > > > > > would someone mind mentioning them?
> > > > > >
> > > > > > Unfortunately, the paper of Choi and Ihm is very detail, and in
> > > > > > the PWCOND code we followed it very closely, so no more
> > > > > > references are needed. In our papers we just extended those ideas
> > > > > > to ultrasoft pseudopotentials, magnetism, spin-orbit coupling and
> > > > > > so on.
> > > > > >
> > > > > > Hope this helped you,
> > > > > > regards, Alexander
> > > > > >
> > > > > > > > Regards,
> > > > > > > > Manoj Srivastava
> > > > > > > > Physics Graduate Student
> > > > > > > > University of Florida,
> > > > > > > > Gainesville,FL, USA
> > > > > > >
> > > > > > > _______________________________________________
> > > > > > > Pw_forum mailing list
> > > > > > > Pw_forum at pwscf.org
> > > > > > > http://www.democritos.it/mailman/listinfo/pw_forum
> > > > > >
> > > > > > _______________________________________________
> > > > > > Pw_forum mailing list
> > > > > > Pw_forum at pwscf.org
> > > > > > http://www.democritos.it/mailman/listinfo/pw_forum
> > > > >
> > > > > Once again thank you very much for your help. I appreciate it.
> > > > >
> > > > > Regards,
> > > > > Manoj Srivastava
> > > > >
> > > > > _______________________________________________
> > > > > Pw_forum mailing list
> > > > > Pw_forum at pwscf.org
> > > > > http://www.democritos.it/mailman/listinfo/pw_forum
> > > >
> > > > _______________________________________________
> > > > Pw_forum mailing list
> > > > Pw_forum at pwscf.org
> > > > http://www.democritos.it/mailman/listinfo/pw_forum
> > >
> > > _______________________________________________
> > > Pw_forum mailing list
> > > Pw_forum at pwscf.org
> > > http://www.democritos.it/mailman/listinfo/pw_forum
> >
> > _______________________________________________
> > Pw_forum mailing list
> > Pw_forum at pwscf.org
> > http://www.democritos.it/mailman/listinfo/pw_forum
>
> _______________________________________________
> Pw_forum mailing list
> Pw_forum at pwscf.org
> http://www.democritos.it/mailman/listinfo/pw_forum