[Pw_forum] Problems with Polarization calculations

[Sohrab Ismail-Beigi]

Dear PWSCF forum,

I've been having a problem running PWSCF for the Berry phase  
calculations of polarization for elements with f (l=3) and higher  
pseudopotentials (i.e. any lanthanides, etc.)  My system is an oxide  
that has La in it, and my pseudopotential is norm conserving and has  
lmax=3.

Here is the error that the program produces upon entering the Berry  
phase code (bp_c_phase.f90):

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      from  ylm_q : error # 181312304
       not programmed for L>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

OK, that is a bit cryptic. The error is coming from the "ylm_q"  
subroutine inside the file "bp_ylm_q.f".  The ylm_q subroutine just  
calculates (real) spherical harmonics and has code for doing so up to  
l=4 (impressive actually).  ylm_q exists with the error above if one  
tries to calculate higher angular momenta than l=4.  The routine is  
given an integer index which it decodes into (l,m):  e.g. index=1  
means l=0; index=4 -> l=1; index = 9 -> l=2, etc.  It is coded up to  
index=25.

Now, the problem is cause by ylm_q being called on being called on  
line 400 of bp_c_phase.f90:

    CALL ylm_q(lmaxq*lmaxq,dk,dkmod,ylm_dk)

It is being asked to compute the angular momentum ylm for a very  
large index lmaxq*lmaxq.   lmaxq itself is set when reading the  
pseudopotential in "allocate_nlpot.f90" where it is set to "lmaxq =  
2*lmaxkb+1".  For lmaxkb=2 (i.e. up to angular momentum l=2), lmaxq=5  
so lmaxq*lmaxq=25 and ylm_q just manages to work.  But for any atom  
with f projectors, lmax=3 -> lmaxq=7, and lmaxq*lmaxq=49 which causes  
the error and failure.

My questions are two-fold:

(1)  First, this is clearly an unintentional bug of some sort.   
Apparently no elements with f projectors (e.g. transition metals)  
were ever used for Berry phase...!??!  Personally, I don't understand  
the "lmaxq*lmaxq" bit at all in this code; it is mysterious.  In  
fact, as far as I can see, the output from ylm_q is being sent to  
"qvan3" which only needs something like lmaxq entries rather  
lmaxq*lmaxq....  But it is hard to be sure without more comments.

(2) On a more global and physics-related point, there are segments of  
the code in bp_c_phase.f90 (in particular the ylm_q calls and related  
quantities) that don't seem to make sense for norm conserving  
pseudopotentials.  Namely, the Berry phase calculation for standard  
orthonormal wavefunctions (from norm-conserving psps) only depends on  
the inner product of Bloch states and neighboring k-points:  one  
shouldn't be needing any angular momentum operators at all!

The ylm_q and related codes, from what I can gather from the code,  
seem to be designed to add a contribution to the electronic part of  
the Berry phase that has the form of some type of non-local atomic  
pseudopotential, which makes me think this is appropriate for ultra- 
soft psps where norm-conservation has been relaxed so one has to have  
some non-local projectors in the augmentation region to fix this up  
when calculating Berry phases.

Therefore, does the bp_c_phase.f90 code implicitly assume ultra-soft  
psps?

Any help you can provide is great --- as it is, there is no way to  
proceed forward unless one finds the real representations of Ylm up  
to some very high l!!!

Thanks,

--Sohrab

---------------------------------------------
Sohrab Ismail-Beigi
Associate Professor of Applied Physics and Physics
Dept. of Applied Physics
Yale University
P.O. Box 208284
New Haven, CT  06520
Phone: 203.432.2107
Fax: 203.432.4283
Email: sohrab.ismail-beigi at yale.edu
http://volga.eng.yale.edu/ - group home page
http://www.eng.yale.edu/aphy/ - dept. page w/ links
http://www.eng.yale.edu/graduate -- info about graduate studies in AP  
and brochures