[Wien] Strange results lapwdm

Lyudmila Dobysheva lyuka17 at mail.ru
Fri Jun 8 09:49:14 CEST 2012


I see in radint.f
            if(krad.le.1)then
           A11(i)=rf1(i,l,is1,if1)
           A21(i)=rf2(i,l,is1,if1)
           A12(i)=rf1(i,l,is2,if2)
           A22(i)=rf2(i,l,is2,if2)
           else if((krad.eq.2).or.(krad.eq.3))then
   ...
           else if(krad.eq.4)then  ! WIEN2k_5 approx.
   ...
           else if(krad.gt.10)then  ! <r**krad>
            A11(i)=rf1(i,l,is1,if1)*rx(i)**(dfloat(krad)/2.)
            A21(i)=rf2(i,l,is1,if1)*rx(i)**(dfloat(krad)/2.)
            A12(i)=rf1(i,l,is2,if2)*rx(i)**(dfloat(krad)/2.)
            A22(i)=rf2(i,l,is2,if2)*rx(i)**(dfloat(krad)/2.)
           else if(krad.lt.-10)then  ! <1/r**krad>
            A11(i)=rf1(i,l,is1,if1)/rx(i)**(dfloat(krad)/2.)
            A21(i)=rf2(i,l,is1,if1)/rx(i)**(dfloat(krad)/2.)
            A12(i)=rf1(i,l,is2,if2)/rx(i)**(dfloat(krad)/2.)
            A22(i)=rf2(i,l,is2,if2)/rx(i)**(dfloat(krad)/2.)
When krad negative it executes first block (krad.le.1), when krad >10 it 
looks like r**(krad), and not 1/r**1

And if to really make the redirection for the case krad < 11 it will do 
1/r**(krad); as krad negative, so again r**(|krad|).
So, I see three points under question
1) in the first block
if(krad.le.1)then
It should be changed to something like:
if(krad.le.1.AND.krad.ge.0)then
2) *rx(i)**(dfloat(krad)/2.) in last but one block to something like:
*rx(i)**(dfloat(krad-10)/2.)
3) /rx(i)**(dfloat(krad)/2.) in the last block to something like:
    /rx(i)**(dfloat( - krad+10)/2.)

It is a very quick view, maybe I didn't understand something.

We need Pavel.

Best regards to all
   Lyudmila Dobysheva

On 07.06.2012 18:58, David Tompsett wrote:
> That is a good idea, with RINDEX= -11 we obtain:
> :XOP  1  0     0.98974     0.00000     0.98974
> On Thu, Jun 7, 2012 at 2:53 PM, Gavin Abo <gsabo at crimson.ua.edu
>     RINDEX = 11 => r**krad for krad > 10 => :XOP 104670.12 (large number)
>     Adding a negative sign will change the calculation:
>     RINDEX = -11 => 1/r**krad for krad < 10 => :XOP ? (small number)
>     File of reference: $WIENROOT/SRC_lapwdm/radint.f
>     On 6/7/2012 7:08 AM, David Tompsett wrote:
>>     Thanks for the response. The option r-index = 11 is not in the UG,
>>     but is documented in Pavel Novak's technical report on the
>>     Hyperfine Field calculation so other users may strike it in future.
>>     On Thu, Jun 7, 2012 at 7:47 AM, Peter Blaha
>>         You have to check directly in the code. Without looking in the
>>         code I do not know what r-index = 11  should do ? I guess this
>>         is not a documented option, but Pavel Novak migth have put in
>>         something ....
>>         Am 06.06.2012 14:11, schrieb David Tompsett:
>>             the radial distribution of orbitals i.e. <r>.
>>             I have attempted the simplest case that I could think of,
>>             the hydrogen atom in a large cell.

>>             ran "x lapwdm -up" with input:
>>
>>             -9.     Emin cutoff energy
>>             1       number of atoms for which density matrix is calculated
>>             1 1 0  index of 1st atom, number of L’s, L1
>>             11 1   r-index, (l,s)-index
>>
>>              From the solution of the Schroedinger equation for a
>>             Hydrogen atom, the
>>             result should be <r> = 1.5 au.
>>             However the result from :XOP in the output of lapwdm is
>>             104670.12. I do
>>             not know what the units are. Also, the value of this
>>             result seems to
>>             vary strongly with the RMT.
------------------------------------------------------------------
Phys.-Techn. Institute of Ural Br. of Russian Ac. of Sci.
426001 Izhevsk, ul.Kirova 132
RUSSIA
------------------------------------------------------------------
Tel.:7(3412) 442118 (home), 218988(office), 250614(Fax)
E-mail: lyu at otf.pti.udm.ru
          lyuka17 at mail.ru
Skype:  lyuka17, lyuka18
http://fti.udm.ru/content/view/25/103/lang,english/
------------------------------------------------------------------


More information about the Wien mailing list